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Preface 


The  purpose  of  this  thesis  was  to  evaluate  the 
performance  of  an  air-to-air  optical  communication  link  in 
the  presence  of  atmospheric  turbulence.  By  varying  certain 
performance  parameters  (i.e.  aircraft  speed  and  altitude),  a 
variety  of  communication  links  can  be  evaluated.  Knowing 
how  a  certain  link  will  perform  can  aid  in  the  determination 
of  the  maximum  reliable  communication  range  between 
aircraft.  I  feel  that  the  results  obtained  from  this 
project  will  provide  some  insight  into  the  limitations  of 
optical  communications  between  aircraft. 

While  writing  this  thesis,  I  find  myself  looking  back 
over  the  months  and  realizing  that  the  contributions  made  by 
others,  whether  great  or  small,  were  an  enormous  help  to  me. 
Recognizing  this,  I  would  like  to  express  my  gratitude  to  my 
advisor,  Dr.  Vaqar  Syed ,  for  his  dedicated  technical 
direction  throughout  the  entire  project.  In  addition,  I 
would  like  to  thank  Majors  Kenneth  G.  Castor  and  Richard  J. 
Cook  for  their  greatly  appreciated  advise  and  counseling. 
Special  thanks  is  due  to  First  Lieutenant  Linden  Mercer, 
AFWAL/AAAI-1 ,  Wr ight-Patterson  AFB,  Ohio,  for  providing  both 
this  thesis  topic  and  an  education  on  the  subject  of 
turbulence  and  optical  equipment  limitations. 

Finally,  I  wish  to  express  my  sincerest  appreciation 
and  gratitude  in  the  dedication  of  this  thesis  to  my  wife, 


Michelle.  Her  support,  love,  and  patience  during  the  period 
of  this  effort  has  made  this  thesis  a  special  success. 
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Jay  Nicholas  Kanavos 


Table  of  Contents 


Page 

Preface .  ii 

List  of  Figures .  vii 

List  of  Tables  .  x 

Abstract  .  xi 

I  Introduction  .  1 

Background  .  1 

Problem  Statement  .  3 

Approach  .  4 

Limitations  .  4 

Organization  .  5 

II  Elements  of  an  Optical 

Communication  System  .  7 

Optical  Communication 

System  Model  .  7 

Physical  Components  .  9 

The  Optical  Transmitter  .  9 

The  Optical  Receiver  .  10 

The  Optical  Channel  .  13 

The  Refractive  Index 

Structure  Constant  .  15 

Probability  Distribution  of 

Intensity  Fluctuations  .  18 

Derivation  of  the  Mean 

and  the  Variance .  21 

Power  Spectral  Density  .  27 

Development  of  the  Denormalized 

Power  Spectral  Density  .  29 

Summary  .  30 

III  Mathematical  Derivation  of 

the  Simulation  .  31 

Introduction  .  . .  31 

Simulation  of  the 

Input  Signal  .  32 

Simulation  of  the 

Atmospheric  Channel  .  33 

Background  Radiation  .  44 

Optical  Beam  Divergence  .  46 

Simulation  of  the  Receiver .  48 


List  of  Figures 


Page 


I 


'  ^ 

r! 

I 


I 


» 


Figure 

1  Block  Diagram  of  an 

Optical  Communication  System  .  . 

2  Receiver  Configuration  . 

2 

3  Experimental  Cn  . 

4  Cn2  vs.  Altitude  . 

5  Communication  System 

block  Diagram  . 

6  Maximum  Propagation  Path  Length  of 

Optical  Beam  Not  in  Saturation  .  . 

7  Lognormal  Distribution  for  a 

Variety  of  Values  . : 

8  Phenomenon  of  Decreasing  Variance 

in  the  Saturation  Region  . 

9  Normalized  Power  Spectral  Density 

of  the  Intensity  Fluctuations  .  . 

10  Denormalized  Power  Spectral  Density 

1 1  Communication  System  Model  .  .  .  . 

12  black  Box  Analogy  . 

13  Block  Diagram  of  the  Gujar 

Kavanaugh  Technique  . 

14  FFT  of  the  Sequences  Z(n) 

and  Y(n)  . 

15  Divergence  of  the  Optical  Beam 

as  the  Path  Length  Increases  .  .  . 

16  Theoretical  Bit  Error  Rate 

for  an  Altitude  of  10,000  Ft  .  .  . 

17  Theoretical  Bit  Error  Rate 

for  an  Altitude  of  20,000  Ft  .  .  . 

18  Theoretical  Bit  Error  Rate 

for  an  Altitude  of  30,000  Ft  .  .  . 


8 

11 

16 

18 

20 

20 

25 

26 

28 

30 

31 

34 

35 

45 

47 

59 

60 

61 


> 


•  *.£yf «■'  ■■vZ-V.  v-  -»•-  -.-  I.  ----  v- 


-V  -3.  -3.  v 


A.  >  •- 


«jM 


vii 


19 


Turbulence  Affected  0] tical  Link 
Altitude:  10,000  ft 

Aircraft  Speed:  425  knts 
Data  Rate:  20,000  bps  .... 

20  Turbulence  Affected  Optical  Link 

Altitude:  10,000  ft 

Aircraft  Speed:  475  knts 
Data  Rate:  20,000  bps  .... 

21  Turbulence  Affected  Optical  Link 

Altitude:  10,000  ft 

Aircraft  Speed:  425  knts 
Data  Rate:  40,000  bps  .  .  .  . 

22  Turbulence  Affected  Optical  Link 

Altitude:  20,000  ft 

Aircraft  Speed:  425  knts 

Data  Rate:  20,000  bps  . 

23  Turbulence  Affected  Optical  Link 

Altitude:  20,000  ft 

Aircraft  Speed:  475  knts 

Data  Rate:  20, 000  bps  . 

24  Turbulence  Affected  Optical  Link 

Altitude:  20,000  ft 

Aircraft  Speed:  425  knts 

Data  Rate:  40,000  bps  . 

25  Turbulence  Affected  Optical  Link 

Altitude:  30,000  ft 

Aircraft  Speed:  425  knts 

Data  Rate:  20,000  bps  . 

26  Turbulence  Affected  Optical  Link 

Altitude:  30,000  ft 

Aircraft  Speed:  475  knts 

Data  Rate:  20,000  bps  . 

27  Turbulence  Affected  Optical  Link 

Altitude:  30,000  ft 

Aircraft  Speed:  425  knts 

Data  Rate:  40,000  bps  . 

28  Representation  of  Possible 

Turbulence  Factor  Sequences  .  . 


64 


65 


66 


67 


68 


69 


70 


71 


72 

75 


viii 


29 


Representation  of  the  Turbule 
Factor  Sequences  for  an  Alti 
of  20,000  ft  and  a  70  km  Pat 
Length  . 


ix 


'  *'  -  '*■ 


Table 


List  of  Tables 


Advantages  and  Disadvantages  of 
an  Unguided  Optical  Communication 
Syst/M . 

Chebychev  Filter  Orders  for  a 
Variety  of  Optical  Propagation 
Environments  . 

Parameters  of  the  Communication 
Links  Used  . 

Transmitter  and  Receiver 
Characteristics  . 

Maximum  Ranges  for  the  Baseline 
Link  and  the  Turbulence  Affected 


Abstract 


This  report  presents  an  analysis  of  the  performance  of 
an  air-to-air  optical  communication  link  in  the  presence  of 
atmospheric  turbulence.  As  aircraft  travel  through  the 
atmosphere,  they  encounter  regions  of  atmospheric 
turbulence.  While  harmless  to  radio  frequency  ( RF ) 
communications,  these  regions  of  turbulence  can  cause  both 
intensity  and  phase  fluctuations  within  an  optical  beam.  As 
a  result,  the  communication  link  bit  error  rate  rises. 

To  evaluate  the  performance  of  such  a  link,  a  computer 
simulation  was  developed.  By  varying  such  parameters  as  the 
speed  of  the  aircraft,  its  altitude,  and  propagation  path 
length,  a  determination  could  be  made  about  link 
performance . 

The  results  obtained  from  the  simulation  showed  that 
atmospheric  turbulence  plays  a  significant  role  in 
determining  link  performance.  It  was  found  that  the  maximum 
reliable  communication  ranges,  of  those  links  affected  by 
turbulence  (real  life  situation)  were  degraded  approximately 
50  -  60  percent  over  those  links  not  affected.  It  was 

also  determined  that  reasonable  changes  in  aircraft  speeds 
had  no  significant  impact  on  the  bit  error  rates,  while  an 
increase  in  altitude  greatly  increased  the  maximum  reliable 
communication  ranges. 


The  Effects  of  Atmospheric  Turbulence  on  an  Air  to  Air 
Optical  Communication  Link 

1.  Introduction 


Background 


The  attention  given  to  optical  communication  systems 


over  the  past  decade  has  been  staggering.  There  has  been  a 
tremendous  effort  under  way  to  develop  workable  system 
designs  for  potential  communication  links.  To  do  this, 
researchers  have  found  that  understanding  the  detrimental 
effects  that  the  atmosphere  has  on  optical  wave  propagation 
is  essential  in  determining  the  "condition"  of  the  received 
wave.  To  date,  researchers  have  attempted  to  describe  the 
effects  which  influence  an  optical  beam  in  both  terrestrial 
and  earth-to-space  satellite  links.  However,  little  work 


has  been  done  in  evaluating  the  effects  which  influence  an 


air-to-air  optical  communication  link. 

As  an  aircraft  travels  through  the  atmosphere,  it 
encounters  regions  of  turbulence.  Turbulence  can  be 


described  as  the  violent  mixing  of  one  substance  wi.  > 
another  substance.  In  this  case,  these  substances  are  the 


ambient  layers  of  the  earth’s  atmosphere.  Throughout  the 
entire  atmosphere,  the  prevailing  winds  constantly  mix 
warmer  air  with  cooler  air.  The  result  is  a  potpourri  of 
thermal  air  pockets,  or  eddies,  flowing  within  the 


atmosphere.  The  temperature  variations  between  these  eddies 
are  very  small  (on  the  order  of  .1-1  degree  C);  however, 
they  are  large  enough  to  cause  index  of  refraction 
fluctuations.  Since  each  eddy  consists  of  a  different 
temperature,  its  index  of  refraction  will  vary.  Thus,  each 
eddy's  index  of  refraction  will  consist  of  a  slightly 
different  value  than  that  of  its  neighbor.  Although  these 
differences  are  typically  small  (on  the  order  of  10“^),  they 
can  have  a  significant  impact  on  the  optical  beam. 

As  an  optical  beam  passes  from  one  eddy  to  another,  the 
changes  in  the  indices  of  refraction  cause  the  beam  to  bend. 
Since  this  beam  passes  through  many  eddies,  parts  of  the 
beam  may  travel  through  more  eddies  than  other  parts  of  the 
beam.  As  a  result,  the  optical  beam  will  arrive  at  the 
detector  with  both  intensity  and  phase  fluctuations. 

These  fluctuations,  known  as  scintillations,  vary 
randomly  with  time.  Thus,  the  turbulence-induced  intensity 
scintillations  become  stochastic  processes.  Since  these 
scintillations  are  not  deterministic,  they  must  be  analyzed 
using  statistical  mathematics.  Therefore,  in  order  to 
determine  the  "condition"  of  the  received  optical  beam,  the 
statistics  of  these  scintillations  must  first  be  determined. 

In  determining  these  statistics,  there  are  three 
quantities  which  must  be  specified:  the  probability 
distribution,  its  mean,  and  its  variance.  It  is  generally 
agreed  that  the  probability  distribution  of  the  intensity 


fluctuations  assumes  a  lognormal  distribution.  However, 
this  assumption  can  only  be  valid  when  the  optical  beam  is 
not  operating  within  the  so  called  "saturation  region."  If 
the  beam  is  operating  within  this  region,  then  its  intensity 
fluctuations  will  no  longer  exhibit  lognormal  statistics. 

In  fact,  researchers  have  yet  to  determine  what  statistics 
are  valid. 

The  lognormal  distribution  has  two  parameters:  its 
mean  and  variance.  It  turns  out  that  the  mean  is  actually 
equal  to  the  negative  of  its  variance.  Thus,  this 
particular  distribution  is  really  a  function  of  only  one 
variable;  its  variance.  The  variance,  however,  is  dependent 
on  a  variety  of  factors  indigenous  to  the  environment 
surrounding  the  transmitted  optical  beam.  Such  factors 
include:  the  propagation  path  distance,  altitude,  and  time 
of  day.  Once  this  variance  is  computed,  the  statistics  of 
the  scintillations  can  be  completely  described. 

Problem  Statement 

Atmospheric  turbulence  induces  both  intensity  and  phase 
fluctuations  in  the  optical  beam.  These  fluctuations  will 
play  a  major  role  in  degrading  communication  linK 
performance,  i.e.  probability  of  bit  error  in  digital 
communication  systems  will  increase.  As  a  result,  mission 
essential  optical  communication  links  could  be  rendered 
inoperative  under  certain  conditions.  To  determine  the 


reliability,  and  feasibility,  of  an  air-to-air  optical 
communication  link,  it  will  be  necessary  to  determine  those 
conditions,  and/or  system  parameters,  which  will  support 
reliable  communications. 

Approach 

This  thesis  will  be  directed  towards  a  simulation  of 
the  atmospheric  effects  which  would  influence  an  air-to-air 
optical  communication  link.  In  determining  how  this 
communication  link  will  be  affected,  the  bit  error  rates 
will  be  determined  at  various  altitudes,  propagation  path 
lengths,  and  aircraft  speeds.  To  do  this,  a  computer  model 
which  simulates  the  transmission  and  detection  of  a  binary 
signal  through  the  turbulent  optical  channel,  will  be 
developed . 

Limitations 

The  primary  thrust  of  this  thesis  will  focus  on  the 
effects  of  scintillation  due  to  clear  air  turbulence.  The 
following  conditions  will  be  assumed  throughout  this 
document : 

1.  Due  to  frequency  stability  factors  inherent  in 
local  oscillators  aboard  aircraft,  a  direct  detection  scheme 
will  be  used.  As  a  result,  intensity  fluctuations  (rather 
than  phase  fluctuations)  of  the  optical  beam  will  be  of 
primary  concern. 
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2.  The  transmitted  optical  signal  will  consist  of 


a  sequence  of  I’s  and  0's  randomly  transmitted  through  the 
atmospheric  channel.  This  will  simulate  a  binary 
information  signal. 

3-  The  optical  beam  will  exhibit  plane  wave 
characteristics .  This  assumption  limits  the  accuracy  of  the 
final  results,  but  greatly  reduces  the  mathematical 
complexity  of  the  problem. 

4.  The  optical  beam  between  the  two  aircraft  will 
be  perpendicular  to  both  aircraft.  Also,  both  aircraft  will 
be  assumed  to  be  at  equivalent  altitudes. 

Organization 

This  thesis  will  be  organized  into  five  chapters  and 
two  appendices. 

Chapter  II  will  include  a  brief  tutorial  on  optical 
communication  systems.  The  characteristics  of  the  turbulent 
channel  will  be  discussed,  as  well  as  a  mathematical 
derivation  of  its  statistics. 

Chapter  III  will  describe  the  approach  taken  in 
producing  the  computer  model,  as  well  as  the  mathematical 
theory  which  justifies  each  of  the  various  steps. 

Chapter  IV  will  present  the  results  obtained  from  the 
simulation.  These  results  will  be  based  on  a  random 
sequence  of  I's  and  0?s  transmitted  through  the  optical 
channel.  An  assessment  of  how  the  system  probability  of 
error  is  affected  by  atmospheric  turbulence  at  various 


altitudes,  propagation  path  lengths,  and  aircraft  speeds 
will  also  be  included. 

Chapter  V  will  consist  of  both  conclusions  and 
recommendations  that  the  author  feels  could  improve  the 
simulation . 

The  appendices  will  consist  of  the  data  points  used  to 
generate  the  performance  curves,  and  the  Fortran  computer 
model  program  listing. 


II.  Elements  of  an  Optical  Communication  System 


.Optical  Communication  System  Model 

The  use  of  optical  communication  techniques  has  long 
been  recognized  as  an  effective  means  of  communicating. 
Optical  communication  techniques  have  advanced  from  the  use 
of  beacon  fires  to  shipboard  blinker  systems  to  the  laser 
communication  systems  of  today.  Present  laser  communication 
systems  have  developed  to  such  a  point  that  they  are  not 
only  practical  for  many  applications ,  but  offer  some 
significant  advantages  over  common  radio  frequency  (RF) 
communications.  However,  optical  communication  links  are 
not  always  the  best  choice  for  a  given  type  of  environment. 
(«T  The  advantages  and  disadvantages  of  unguided  optical 

communication  systems  are  given  in  Table  I. 

The  block  diagram  of  a  generalized  optical 
communication  system  model  is  shown  in  Figure  1.  In  this 
model,  the  information  signal  can  be  thought  of  as  a 
sequence  of  discrete  symbols  acquired  from  either  a  sampled 
waveform  or  a  digital  source.  This  discrete  input  signal  is 
then  modulated  onto  an  optical  carrier.  The  carrier  is  then 
transmitted  as  an  optical  light  field,  or  beam,  through  the 
turbulent  optical  channel.  At  the  receiver,  the  field  is 
optically  collected  and  processed,  along  with  any 
interference  and  noise.  Thus,  the  output  signal  is  a 
_  corrupted  version  of  the  input  signal. 
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fable  I 

Advantages  and  disadvantages  of  an  unguided  optical 
communication  system  (1:34) 


Advantages 


High  data  rate  capability 


Provides  a  good  level  of 
security 


Low  power  requirements 


Relatively  immune  to 
crosstalk 


Small  antenna  size 


No  communication  license 
required 


Exploits  unused  part  of 
electromagnetic  spectrum 


INPUT 

SIGNAL 


Mnnm a tor 

Disadvantages 


Uncertain  reliability  due 
to  atmospheric  conditions 


Requires  very  accurate 
pointing  and  tracking 


Less  suited  to  broadcasting 
because  of  narrow  beam 


Slightly  higher  level  of 
noise  in  the  received 
signal  due  to  quantum 
nature  of  detection 


Possible  eye  hazard 


OPTICAL 

CHANNEL" 


DETECTOR  H  DEMODULATOR 


OUTPUT 

SIGNAL 


Figure  1.  Block  diagram  of  an  optical  communication 
system. 
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Since  this  thesis  is  concerned  with  the  effects  of 
atmospheric  turbulence,  the  bulk  of  the  discussion  in  this 
chapter  will  concentrate  on  the  optical  channel.  However, 
to  understand  the  entire  optical  communication  system  model, 
a  basic  understanding  of  its  physical  components  is  needed. 

Physical  Components 

The  main  equipment  components  in  the  optical 
communication  system  are  the  transmitter  and  the  receiver. 


The  Optical  Transmitter 

The  optical  transmitter  can  be  represented  as  a  cascade 
connection  of  an  optical  source,  modulator,  and  antenna 
(lens  system). 

Optical  sources  (such  as  lasers  and  light  emitting 
diodes  (LED's))  produce  the  optical  electromagnetic  fields 
which  act  as  the  carrier  for  t-he  information  signal.  Among 
some  of  the  most  popular  laser  sources  used  today  are: 

Nd:YAG  (1.06  Mm),  Nd:YAG  (.53 Mm),  and  CC^  (10.6Mm  and 
9.6  Mm ) • 

In  order  to  efficiently  transmit  information  from  one 
location  to  another,  the  information  signal  must  first  be 
conditioned  in  some  manner.  This  conditioning  is  most 
commonly  referred  to  as  modulation.  Modulation  can  be 
described  as  the  process  of  combining  an  information  signal, 
with  a  carrier,  to  construct  a  more  suitable  signal  form  for 
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transmission.  There  are  a  variety  of  optical  modulation 
techniques  available  today.  Among  some  of  the  schemes  used 
are:  intensity  modulation  (IM),  phase  modulation  (PM),  and 
polarization  modulation  (PL).  In  choosing  which  modulation 
scheme  is  best  for  a  given  type  of  system,  there  are  several 
factors  which  must  first  be  considered.  For  example,  the 
kind  of  data  to  be  transmitted  and  the  type  of  detection 
equipment  to  be  used  are  important  considerations.  If  a 
direct  detection  receiver  scheme  were  to  be  used,  then  some 
type  of  intensity  modulation  would  be  a  logical  choice. 

Once  the  system  "environment1'  has  been  defined,  then  the 
choice  of  a  modulation  scheme  can  be  made. 

The  Optical  Receiver 

The  purpose  of  the  receiver  is  to  collect  the  incoming 
field  of  the  optical  beam.  This  is  done  by  using  a  series 
of  lenses  which  focuses  the  incoming  field  onto  a 
photodetecting  surface. 

As  can  be  seen  in  Figure  2,  the  receiving  lens  is 
located  in  the  aperture  plane,  which  has  an  area  "Ar"  and  a 
focal  length  "f."  The  lens  focuses  the  received  optical 
field  onto  the  focal  plane  behind  the  lens.  The  focused 
field  appears  as  a  diffracted  pattern  in  the  focal  plane. 

The  photodetector  is  located  in  the  focal  plane,  and 
responds  to  the  portion  of  the  diffracted  field  on  its 
surface.  The  photodetector  converts  the  field  intensity  of 


Figure  2.  Receiver  configuration  (2:27). 

the  optical  carrier  into  an  electrical  signal  (2:27). 

Typically,  photodetectors  fall  into  two  classes: 
thermal  and  photon  detectors.  Since  thermal  detectors 
inherently  have  low  bandwidth  capabilities,  they  will  be  of 
little  interest  to  us.  Photon  detectors,  however,  offer 
significant  advantages  over  thermal  detectors,  and  can  be 
subdivided  into  classes  according  to  their  method  of 
optical-to-electrical  conversion.  The  most  common 
conversions  are:  photoemission  effect,  photoconductive 
effect,  photovoltaic  effect,  and  the  photoelectromagnetic 
effect . 

The  photoemissi ve  effect  involves  the  emissions  of 
electrons  from  a  vacuum  tube  cathode  in  response  to  optical 
excitation,  while  the  latter  three  types  of  photo-effects 


are  all  associated  with  semiconductor  detectors.  Here,  the 
absorption  of  photons  on  the  detector  surface  leads  to  a 
change  in  the  concentration  of  charge  carriers  in  the 
material.  This  change  in  the  concentration  of  charge 
carriers  produces  the  signal  current. 

After  the  photodetector  has  converted  the  optical  beam 
into  an  electrical  current,  this  current  is  then  processed 
by  some  type  of  signal  processor.  Unfortunately,  this 
processed  signal  is  not  the  same  signal  that  was 
transmitted.  Disregarding  the  fact  that  the  signal 
possesses  turbulence-induced  intensity  fluctuations,  it  has 
also  been  corrupted  by  both  thermal  and  shot  noise,  as  well 
as  dark  current. 

Dark  current,  thermal,  and  shot  noises  are  all 
corruptive  signals  which  originate  within  the  receiver 
components.  Each  is  random  by  nature  and  has  its  own  set  of 
statistics.  All  three  forms  of  these  noises  are  present  in 
the  photodetector,  however,  there  are  techniques  which  can 
be  used  to  lessen  their  effects  on  a  signal.  For  example, 
since  the  dark  current's  power  spectra  consists  of  an 
impulse  located  at  f=0,  it  can  be  easily  removed  by  a  DC 
block  filter. 

Unfortunately,  the  effects  due  to  thermal  and  shot 
noise  cannot  be  eliminated  as  easily.  The  best  that  can  be 
done  is  to  limit  the  effects  of  one  or  the  other:  i.e.  allow 
one  to  dominate  over  the  other.  Doing  this  allows  the 


potential  system  designer  to  use  the  statistical  quantities 
of  the  dominating  noise  in  determining  the  system 
probability  of  error  and  signal  to  noise  ratios.  This 
greatly  simplifies  the  calculations,  yet  still  provides  a 
reasonable  amount  of  accuracy  in  the  determination  of  these 
quantities . 

The  Optical  Channel 

Temperature  instability,  within  the  earth's  atmosphere, 
(resulting  from  rising  warmer  air  mixing  with  descending 
cooler  air)  produces  a  state  of  turbulence.  This  turbulence 
causes  fluctuations  in  the  index  of  refraction  within 
various  neighboring  eddies.  Equation  (1)  (3:10)  is  the 
standard  relation  used  to  determine  the  value  of  the  index 
of  refraction  at  a  particular  point  in  both  time  and  space. 

n  =  77.6(1+(7.52x10"3)(  A ~2) ) ( P/T) (1 0"6)+1  (1) 

where 

n  =  index  of  refraction 
A  =  wavelength  of  light 
P  =  pressure  (millibars) 

T  =  temperature  (°  K) 

By  observing  Eq .  (1),  it  can  be  seen  that  as  either  the 
pressure  or  the  temperature  changes,  so  too  does  the  index 
of  refraction.  As  a  result,  eddies  of  varying  indices  of 
refraction  are  formed.  Even  though  the  sizes  of  these 
eddies  can  vary,  researchers  have  found  it  convenient  to 
classify  them  into  two  groups.  These  groups  are  known  as 
outer  and  inner  scale  variations. 


The  outer  scale  variations,  denoted  L  ,  is  the  name 

o 

used  to  represent  those  groups  of  eddies  which  are 
physically  larger  in  size.  Values  of  Lq  depend  mainly  on 
the  height  of  these  eddies  above  the  ground.  Typical  values 
for  Lq  are  100  meters  or  1/5  the  height  above  the  ground  - 
whichever  is  less  (4:1527)- 

If  there  is  a  great  deal  of  turbulence  present  in  the 
atmosphere,  then  some  of  the  outer  scale  variations  can 
become  unstable  and  begin  to  break  apart.  As  a  result,  they 
become  smaller.  The  resulting  smaller  eddies  are  called 
inner  scale  variations.  Inner  scale  variations  are  usually 
only  a  few  millimeters  in  size  and  denoted  as  1  (4:1526- 

1527).  To  get  a  feel  for  the  size  of  these  inner  scale 
variations,  their  dimensions  can  be  approximated  by  the 
following  relation  (5:29): 

10*  L0/(h3/,) 

where 

R  =  Reynold's  number 

Since  the  atmosphere  is  dynamic,  the  sizes  and  number 
of  these  outer  and  inner  scale  variations  are  constantly 
changing.  Thus,  as  the  optical  beam  propagates  through  the 
atmosphere,  the  changes  in  the  indices  of  refraction  cause 
scintillations  within  the  beam.  To  determine  the 
"condition''  of  this  received  optical  beam,  these 
scintillations  need  to  be  predicted.  Since  they  vary  with 
time,  statistical  mathematical  techniques  must  be  used. 
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To  describe  the  scintillations  statistically,  there  are 
three  parameters  which  must  first  be  determined.  They  are: 


.the  probability  distribution  function  of  the  intensity 
fluctuations,  its  mean,  and  its  variance.  These  three 
quantities,  as  well  as  the  refractive  index  structure 
constant  and  the  power  spectral  density,  will  all  need  to  be 
determined  in  order  to  predict  the  effects  of  turbulence  on 
the  optical  beam. 

The  following  sections  will  provide  the  necessary 
mathematical  derivations  for  the  determination  of  the  afore 
mentioned  quantities. 


The  Refractive  Index  Structure  Constant 

The  refractive  index  structure  constant,  denoted  as 

2 

Cn  ,  is  a  measure  of  how  strong  the  index  of  refraction 

fluctuations  are.  Since  these  fluctuations  are  a  function 

of  temperature  variations,  then  it  would  seem  reasonable 
2 

that  Cn  must  also  vary  with  the  temperature.  In  fact,  it 
does ! 

2 

Measurements  (3:20-23)  of  Cn  have  been  made  by  using 

airborne  mounted  or  balloon  borne  sensors  which  measure  the 

temperature  at  various  points  in  the  atmosphere.  This 

temperature  data  is  then  used  to  compute  a  temperature 

2 

structure  parameter,  Cfc  which  is  given  (3:12-13)  by 


.-2/3 


where 


E{  [T(r1)-T(r2)]£  } 


r  =  distance  between  r.  and  r~ 
T(r^),T(r2)  =  temperature  at  either  r^  or  r. 


'  '  '.i 

- 

9  i 
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Using  the  results  from  Eq .  (2),  Cn  can  then  be  calculated 
(3:13) 

Cn2  =  [ 7 9 x IQ"6  (P/T2)]2  Ct2  (3 

where 

P  =  pressure 
T  =  temperature 

By  using  Lq .  (3),  researchers  could  take  experimental 

2 

temperature  data  and  compute  numerical  values  for  Cn  at  a 

2 

particular  point  in  time.  Data  obtained  from  averaged  Cn 
measurements  can  be  seen  in  Figure  3.  The  "hump"  at  40,000 
ft,  in  the  figure,  represents  values  of  Cn  in  the 
tropopause . 

p 

Using  Figure  3,  approximate  values  of  Cn  may  be 
calculated.  By  taking  altitude  intervals,  and  then  fitting 


I  ■ 


I 


curves  to  the  Cn  data,  the  following  expressions  can  be 
derived  for  determining  Cn  . 


10  m  <  altitude  <  1  km  : 


Cn2  =  [ [ 5 . 3x 1 0“8  -  5 


.  07x10"8/altitude-1o\]  (4.642)1 
\  990  )]  1 


1  km  <  altitude  <  8  km 


Cn2  =  [  [2 .  3x  1  0“9  -  1 . 97x10~9  ^altitude- 1  000  (4.642)j 


8  km  <  altitude  <  1 0  km  : 


Cn2  =  |  1 3 - 3x 1 0  — 1 0  -  3x10_1 1  /altitude-80001 
1  1  ^  1000  / 

+  3x10-11  /altitude-8000)  / altitude-8000 

l  TMO  7  y  TXRH5  ~ 


-  1^](  4. 642)] 


10  km  <  altitude  <  14  km  : 


If3.3x10~10  +  3-  17x10"9  /altitude-IOOOC) 

1 1  l  2000  / 

3.  17x10“^  /altitude-10000) /altitude-10000 

y  zmg  )  y 


-  lj](4.642)] 


Using  the  expressions  above,  different  values  for  Cn^  can  be 

found  by  substituting  the  appropriate  altitude  into  the 

2 

appropriate  expression.  Figure  4  illustrates  Cn  for  a 
variety  of  altitudes. 


.  %  «**  • 
'/.•'aV-'*.* 
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Figure  4.  Cn  vs.  altitude 

Probability  Distribution  Function  of  the  Intensity 
Fluctuations 

The  effects  of  atmospheric  turbulence  on  an  optical 
wave  are  unlike  those  of  the  corruptive  noise.  The  effects 
due  to  turbulence  are  usually  represented  as  multiplicative 
(7:1628);  while  those  effects  caused  by  noise  are  additive. 
A  block  diagram  representing  this  convention  is  shown  in 


Figure  5.  In  communications  theory,  these  multiplicative 
effects  are  more  commonly  referred  to  as  time  varying 
fading.  Thus,  as  an  optical  beam  propagates  through  the 
atmosphere,  its  intensity  will  fluctuate,  or  fade,  with 
time . 

Experimentation  has  shown  that  the  intensity 
fluctuations  of  the  optical  beam  appear  to  exhibit  lognormal 
probability  statistics.  However,  this  is  only  true  if  the 
optical  beam  is  not  operating  within  the  saturation  region. 

A  horizontal  propagating  beam  is  said  to  be  in  the 
saturation  region  when  the  following  condition  is  satisfied: 


k7/6  l11/6  Cn2  >  5 


(A) 


where 


k 

L 


wavenumber 

path  length  of  the  beam 


Using  this  relation,  Figure  6  illustrates  the  maximum 
propagation  path  length  as  a  function  of  altitude.  As  long 
as  this  path  length  is  below  the  curve  ,  in  Figure  6,  then 
the  optical  beam  is  not  considered  to  be  in  the  saturation 
region . 

While  in  the  saturation  region,  there  has  been  some 
doubt  as  to  whether  the  statistics  of  the  intensity 
fluctuations  remain  lognormal.  In  fact,  there  has  been  a 
great  deal  of  research  conducted  which  has  attempted  to 
determine  what  type  of  statistics  are  valid.  Among  some  of 
the  statistical  distributions  that  have  been  considered  are 
the  Rayleigh,  Rice-Nakagami ,  and  K  distributions. 
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Figure  5-  Communication  system  block  diagram 
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Figure  6.  Maximum  propagation  path  length  of 
optical  beam  not  in  saturation 


Since  it  has  been  shown  that  the  probability 
distribution  of  the  intensity  fluctuations  exhibits 
lognormal  statistics,  then  all  that  needs  to  be  determined 
are  the  parameters  of  this  distribution,  which  are,  its 
v^’iance  and  its  mean. 


Derivation  of  the  Variance  and  the  Mean 

The  variance  and  the  mean  of  the  lognormal  probability 
distribution  can  be  found  from  using  statistical  mathematics 
in  the  following  manner  (8). 

From  elementary  statistics,  it  is  known  that  if  a 
normal  random  variable  (  x  )  is  written  as 

A  =  exp(  X  ) 

then  A  is  said  to  be  lognormally  distributed.  Since  the 
received  field  of  the  optical  beam  is  distributed 


lognormally,  then  it  can  be  represented  as 


where 


U  =  UQ  exp (  X  +j  6 ) 


field  of  the  received  beam 
field  of  the  transmitted  beam 
log  amplitude 
phase 


Then  the  intensity,  I,  of  the  received  field  is  given  by 

I  =  I  U  1 2 

=  IUQ|2  exp ( X  +j  6 )exp( x  -j  0  ) 
=  IQ  exp (2  x  )  i 


I 


Or,  I/I0  =  exp ( 2 x  )  (6) 

If  we  take  the  natural  log  of  both  sides  of  Eq .  (6),  we  get 

In ( I/IQ )  =  ln(exp(2  X  )) 

=  2X  (7) 

But  In ( I/I  )  =  ln(A2/A  2) 

o  o 


=  21n ( A/A  ) 
o 


where 


A  =  amplitude  of  the  received  wave 
Aq  =  amplitude  of  the  transmitted  wave 


Combining  Eqs.  (7)  and  (8)  we  obtain 

21n ( A/A  )  =  In ( I/I  ) 
o  o 

=  2  X 

Or,  ln(A/AQ)  =  x  =  log  amplitude 


Since  the  normalized  intensity,  I/1Q,  is  lognormally 
distributed,  then  the  probability  distribution  function, 
p(l/IQ)  is  given  by 


P(I/Io)  =  1 


(I/I  0)/F?o 


ln(I/Io)  -  E {In ( I/IQ  )  J ( 9  ) 


To  determine  O  ,  observe  that  it  is  defined  by 

a2  =  E{ln(I/lo)2)-E2{ln(I/lo)} 


Using  Eq.  (7) 


a2  =  E{4  X2}-E2{2  X  } 

:  U{  X  2  }  -4  E2{  X  } 
=  4  [E{  X 2 } — E 2 {  X  }] 


»*•  .jiili .  *  ■*’  «*■  .*•  -  .*•  #'•  /•  .‘•V-V  >  .*•  *-  .*• 


(10) 


a 


2 


Thus,  the  variance  of  the  log  intensity  is  equal  to  H  times 
the  variance  of  the  log  amplitude. 

To  determine  the  mean  or  E{ln(l/I  )}  in  Eq.  (9),  the 
physical  Law  of  Conservation  of  Energy  must  be  used.  In 
short,  this  means  that  the  average  amount  of  energy  crossing 
the  transmitter  plane  has  to  equal  the  total  amount  of 
energy  crossing  the  receiver  plane. 

Thus 

E{|Urx|2}  =  £{|Utx|2} 

Or,  E{Irx}  =  E{Itx}  (11) 


where 


|U  |_  =  I  =  intensity  at  the  receiver 
*Utx*^  =  tx  =  intensity  at  the  transmitter 

Taking  the  expected  value  of  Eq .  (5) 

E { I }  =  E{Io)  E {exp (2 X  )} 


Dividing  both  sides  by  E { IQ } ,  and  using  Eq .  (11) 


1  =  E {exp( 2  X  ) } 


E { exp ( 2  X )} 


*  / 

•'-CD 


exp(2  x  )  p( x  )  dx 


Since  x  is  normally  distributed: 


E{exp(2  x  ) } 


-f 


exp(2  x  )  1  exp  -(  x-  X  )2  dx 


23 


1  exp 


Using  the  relation 


exp(-ax  +  bx)  dx  = 


VT  exp(4) 


The  last  equation  above  reduces  to 


E  {exp( 2  X ) }  = 


=  exp  - 


1  exp/ -X2NU/2T  exp  fl2  +  g?) 

^  v^r' 

p/-x2  Wpf/  <Tx\/4  +  4  x  +  x2\ 
\2af)  [V  2  A  HI  n?)\ 


=  exp (2  +  2  X  ) 


where 

X=  E  {  X  } 

Thus  E{exp(2X)}  =  exp(2  Oy  +  2X) 

Now 

E  {exp  ( 2  X  )]  =  exp  (2  ff2  +  2X  )  =  1 

Thus 

*  =  -«x2  <12 

E  {In  ( I/IQ) )  =  E{2X  } 

=  2  E{  X  } 

=  -2  O*  (13 


Therefore,  the  revised  probability  distribution  function  of 
the  intensity  fluctuations  can  now  be  determined  by 


Figure  7.  Lognormal  distribution  for  a  variety  of  a. 
values. 

substituting  Lqs.  (10),  (12),  and  (13)  into  Eq  .  (9)  to  get 

pd/I  )  a  1  exp/-(  fcln(l/Io)  ♦  ol  )2\  (14) 

(T7TjJT*\  \  ) 

Figure  7  shows  a  plot  of  the  probability  density  function 

2 

p(l/I  )  as  a  function  of  O .  . 

O  A 

2 

As  the  inequality  in  Eq .  (4)  gets  closer  to  5,  <?x 

increases  to  a  value  of  approximately  .64  .  However,  when 

2 

the  beam  "saturates"  (condition  in  Eq .  (4)  is  met),  then 
no  longer  increases.  In  fact,  it  has  been  noted  that  as 

p 

either  the  path  length  or  Cn  increases,  there  is  a  slight 
decrease  in  .  This  phenomenon  is  illustrated  in 
Figure  8. 


saturation 


100's  -  1000's  of  meters 


Figure  8.  Phenomenon  of  decreasing  variance 
in  the  saturation  region. 


Since  the  probability  distribution  function  of  Eq .  (14) 

2 

is  a  function  of  only  one  variable-  <TX  ,  then  this  parameter 
must  be  determined.  In  his  works,  Tatarski  has  shown  that 
the  variance  of  the  log  amplitude  can  be  expressed  (4:1534) 


at  =  2w2  k2  L 


!}- 


k  sin  /  T  L 
TT  V~k 


ih)] 


$ (T)  T  dT 


where 


k  =  wavenumber 

L  =  path  length 

T  =  spatial  wavenumber  (2T/1) 

1  =  scale  size  of  an  eddy 

4>(T)=  power  spectral  density  of  the 

index  of  refraction  fluctuations 


However,  Tatarski  has  also  shown  that  if 


l  <<  JTl 


where 


size  of  an  inner  scale  variation 


(17) 


Then  Eq .  (15)  can  be  approximated  by 


124  Cn2  k7/6  L11/6 


Since  Eq .  (16)  is  always  true  for  reasonable  path  lengths 
(lengths  greater  than  10  meters),  then  Eq .  (17)  can  be  used 
with  confidence. 


Power  Spectral  Density 

The  frequency  spectrum  of  the  intensity  fluctuations  is 
known  to  be  caused  by  the  motion  of  the  atmosphere.  This 
motion  can  be  separated  into  three  groups: 

1)  Mean  motion  of  the  wind 

2)  Changes  in  the  wind  direction 

3)  Internal  "mixing"  motions  of  the  atmosphere 
due  to  turbulence 

Since  this  thesis  is  concerned  with  communications  between 
two  aircraft,  then  it  will  be  assumed  that  the  aircraft  are 
stationary,  and  the  atmosphere  is  moving  rather  than  vice 
versa.  If  the  "Frozen  Turbulence  Assumption"  is  used,  then 
groups  2  and  3>  above,  may  be  neglected.  This  is  due  to  the 
fact  that  the  frozen  turbulence  assumption  assumes  that  any 
variations  in  the  indices  of  refraction,  measured  at  a 
point,  will  be  caused  by  the  uniform  motion  of  the 
atmosphere  past  that  point.  Thus,  internal  motion  within 
the  atmosphere  may  be  neglected  (4:1537).  As  a  result, 
group  1  wind  motions  will  be  of  primary  interest.  This 


frozen  turbulence  assumption  may  only  be  used 

when  vXL  <<  L  ,  which  is  true  most  of  the  time. 

’  o 

Tatarski  has  rioted  that  the  "normalized"  power  spectral 
density  C  U  ( f ) )  of  the  intensity  fluctuations  can  be  written 
(5:218)  as: 

r®  -U. 

U(f)  =  f  W(f)  =  1.35ft  /  [l  -  sinft2  +  ft2)  1<  t2  +  ft2)  at  (18) 
Q  a7  Jl  U2  +  ft2)  J 


where 


W(f)  =  true  power  spectral  density 

ft  =  r/f0 _ 

f0  =  ws/^/2  r  XL 

ws  =  wind  speed 


Typically,  most  literature  represents  the  normalized  power 
spectral  density  as  shown  in  Figure  9.  However,  by 
rewriting  the  axes,  a  "denormalized"  version  of  Figure  9  can 
be  found. 


f0W(f )/  o£ 


-0.5 


In ( f / f J 


Figure  9.  Normalized  power  spectral  density 

of  the  intensity  fluctuations  (5:218). 


Development  of  the  Denormalized  Power  Spectral  Density 


Before  the  denormalized  version  of  the  power  spectral 
density  of  Figure  9  can  be  determined,  values  for 
fQW(  f  )/<jr^  vs .  the  corresponding  values  for  log(f/fQ)  must  be 
found.  This  can  be  done  by  directly  reading  them  from 

Figure  9.  For  example,  when  f«  W  ( f )  =  .43,  then 

Ox 

1 o  g ( f / f o )  =  0  (19) 

Thus , 


f 


f  =  ws 

0  72-FTT 


(20) 


bince  Figure  9  was  based  on  ws  =  Im/s,  L  =  1000m, 
and  X  =  5121  A  then 
f  =  f 


1 


y 2  ir  (5121x10‘10)(1x103) 


r  =  17.63  Hz 


Finally,  to  determine  that  value  of  W(f)  (denormalized  power 
spectral  density  value)  which  corresponds  to  f  = 

9  2 

17.63  Hz,  needs  to  be  calculated.  To  do  this,  Cn  must 
first  be  determined  from  the  available  parameters-  namely 

2 

the  altitude.  By  using  the  relations  derived  earlier  ,  Cn 

2  2 

can  be  calculated.  Once  Cn  has  been  determined,  <?x  is 
found  from  Eq .  (17). 

From  here  on,  W(f)  can  be  found  for  any  value  of 
frequency  in  the  power  spectral  density  by  reading  Figure  9 
and  transforming  them  as  shown  above.  An  example  of  a 
denormalized  power  spectral  density  is  shown  in  Figure  10. 

As  can  be  seen,  Figure  10  actually  resembles  a  Chebychev  low 


pass  filter.  As  will  be  seen  in  Chapter  III,  this  fact  will 
be  useful  in  modeling  the  power  spectral  density. 

Summary 

In  summary,  it  has  been  determined  that  the  intensity 
fluctuations  of  the  optical  beam  exhibit  lognormal 
statistics,  and  are  a  function  of  many  parameters.  To  model 
these  fluctuations,  the  following  expressions  have  been 
derived 

E  {  X  }  =  -  a £ 

0 £  =  .124  Cn2  k7/6  L1 1/6 

pd/I  )  =  1  exp  (  -(  k  lnd/Io)  ♦  cl  )2\ 

(I/I0)/8?<fy  \  2  0^  / 

In  addition,  a  means  of  representing  the  power  spectral 

density  of  the  intensity  fluctuations  has  been  developed, 

which  in  turn  will  be  used  as  a  basis  for  the  computer 

simulation . 


W(f) 


Figure  10.  Denormalized  power  spectral  density. 


III.  Mathematical  Derivation  of  the  Simulation 


Introduction 

In  Chapter  II,  the  intensity  fluctuations  of  the 
optical  wave  were  found  to  be  distributed  lognormally.  It 
was  also  determined  that  these  scintillations  affected  the 
beam  in  a  multiplicative  manner  rather  than  in  an  additive 
one.  To  represent  this  convention,  the  model  shown  in 
Figure  11  will  be  used  as  a  basis  for  this  thesis.  To 
determine  how  reliable  this  communication  system  will  be, 
each  of  the  factors  shown  in  Figure  11  must  be  modeled 
mathematically . 

The  ultimate  goal  of  this  simulation  will  be  to 
transmit  a  binary  digital  signal,  s(n),  through  the 
atmospheric  channel  and  into  the  receiver.  While  in  the 
receiver,  the  received  signal  sequence,  r(n), 
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Figure  11.  Communication  system  model. 
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will  be  corrupted  by  noise  and  then  processed.  The  receiver 
will  use  a  threshold  detection  technique  to  determine 
whether  a  1  or  a  0  was  transmitted.  By  comparing  what  was 
sent  to  what  the  receiver  decided  was  sent,  a  bit  error  rate 
can  be  established.  Running  this  simulation  with  varying 
environmental  conditions  (i.e.  altitude,  propagation  path 
length,  and  aircraft  speed)  should  give  a  good  indication  of 
the  effects  of  turbulence  and  noise  on  an  air-to-air  optical 
communication  link. 

The  ensuing  sections  of  this  chapter  will  present  the 
method  and  approach  used  to  model  each  of  the  elements  shown 
in  Figure  11.  Once  modeled,  the  complete  simulation  can  be 
run  and  the  bit  error  rate  determined. 

Simulation  of  the  Input  Signal 

Since  we  are  dealing  with  a  digital  communication 
scheme  (on/off  keying),  the  input  signal  can  be  represented 
as  a  random  sequence  of  I's  and  0's  being  transmitted 
through  the  optical  channel. 

Using  the  International  Mathematical  and  Statistical 
Library  (IMSL)  subroutines,  pseudorandom  numbers  between  0 
and  1,  can  be  generated.  Once  generated,  these  numbers  can 
be  transformed  into  binary  I’s  and  0's  using  the  following 
logic  scheme 

Let:  x  =  IMSL  pseudorandom  number 


y  =  binary  signal  bit 


If: 

x  >  .  5 

then 

y  =  i 

x  <  .  5 

then 

O 

II 

>» 

By  repeatedly  doing  this 

( about 

750,000 

times  per  simulation 

run),  an  information  signal  sequence  containing  750,000  I's  and 
0's  can  be  transmitted  through  the  optical  channel. 

Simulation  of  the  Atmospheric  Channel 

From  the  communication  block  diagram,  the  optical  channel  is 
shown  to  consist  of  two  elements:  atmospheric  turbulence,  and 
bacKground  noise.  Each  of  these  elements  are  random  processes 
and,  as  such,  must  be  dealt  with  using  statistical  mathematics. 

Using  Figure  1 1  as  a  reference,  the  modeling  goal  for  the 
optical  channel  will  be  to  develop  a  sequence  of  numbers,  t(n), 
which  can  be  multiplied  by  the  signal  sequence,  s(n).  The 
multiplicative  product,  S-jCn),  will  represent  the  corrupted 
signal  sequence.  The  intensity  of  each  bit  in  s^n)  must 
fluctuate  lognormally.  Once  this  sequence  is  produced,  it  can 
then  be  corrupted  by  noise  from  any  background  radiation  present 
at  the  time. 

As  mentioned  earlier,  various  researchers  have  shown  that 
based  upon  theoretical  and  experimental  considerations  the 
intensity  fluctuations  of  the  optical  signal  exhibit  lognormal 
statistics.  However,  when  these  researchers  have  turned  to 
experimental  methods,  the  optical  beams  used  have  always 
consisted  of  a  constant  intensity.  To  grasp  the  approach  taken 
in  this  thesis  for  modeling  the  scintillation  effects,  an  analogy 
between  a  real  experiment  and  an 


Figure  12.  Black  box  analogy. 

illustrative  example  can  be  drawn. 

This  analogy  is  illustrated  in  Figure  12.  Here  a  "digital” 

laser  beam  of  constant  intensity  is  passed  through  a  "black  box." 

Let  the  term  "constant  intensity"  mean  a  continuous  stream  of 

digital  I's.  As  this  beam  passes  through  the  box,  its  intensity 

gets  "multiplied"  by  some  other  sequence  of  numbers  (call  them 

turbulence  factors).  Assume  that  these  turbulence  factors 

» 

fluctuate  lognormally  and  have  a  power  spectral  density 
resembling  the  frequency  response  of  a  low  pass  filter.  As  the 
beams  exits,  its  intensity  is  no  longer  constant.  Instead,  it 
fluctuates  lognormally  and  has  a  low  pass  filter-like  power 
spectral  density.  Since  its  intensity,  entering  the  box,  was 
constant  (i.e.  a  sequence  of  digital  I's),  then  the  sequence  of 


3^ 


bits  coming  out  of  the  box  must  be  equal  to  the 
multiplicative  turbulence  factor  sequence  of  bits  in  the 
box . 

Therefore,  we  can  model  a  turbulence  influenced  optical 
beam  using  the  black  box  analogy.  This  means  that  the 
turbulence  element,  t(n)  in  Figure  11,  can  be  modeled  as  a 
random  sequence  of  numbers  whose  magnitude  fluctuates 
lognormally  and  power  spectral  density  (PSD)  resembles  a 
Chebychev-like  low  pass  filter.  This  sequence  of  numbers 
will  then  be  multiplied  by  the  signal  s(n).  To  generate  a 
random  series  of  numbers  to  meet  these  requirements ,  a 
technique  developed  by  Gujar  and  Kavanaugh  (9:716-719)  can 
be  used. 

The  principle  behind  this  technique  is  illustrated  in 
Figure  13.  A  random  signal,  x(t),  has  a  gaussian 
distribution,  p  (x),  and  a  known  PSD,  S  ( u>  ) .  The  transfer 


Figure  13.  Block  diagram  of  the  Gujar  Kavanaugh 
technique  (10:  17) • 


function,  H(w),  of  the  linear  filter  is  chosen  such  that  it 
will  produce  an  output  signal,  z(t),  which  is  also  gaussian 
but  has  the  desired  PSD.  The  signal,  z(t),  is  then  passed 
through  a  non-linear  device.  The  resulting  signal,  y(t), 
will  have  both  the  desired  probability  distribution  and  PSD. 

The  desired  transfer  function,  H(a>),  can  be  determined 
from  elementary  stochastic  theory  in  the  following  manner 
If:  S^ ( w  )  =  PSD  of  the  input  signal 

Sz(W)  =  PSD  of  the  output  signal 
H (  u>  )  =  filter  transfer  function 


Then  : 

Sz< « )  =  |H( w ) |2  Sx 

If: 

S  ( u»  )  is  white 

X 

Then  : 

Sx(u>)  =  1 

Sz(  u>  )  =  | H(  W  )  |  2 

H(  U)  )  a  ^Sz(  W  ) 

This  means  that  the  magnitude  squared  of  the  filter's 
transfer  function  is  equal  to  the  desired  PSD. 

To  duplicate  the  above  technique  for  the  computer 
simulation,  the  input  signal,  x(t),  must  have  a  gaussian 
distribution  and  a  white  PSD.  This  input  signal  can  be 
generated  from  the  IMSL  computer  subroutines.  Since  the 
gaussian  random  variables  are  generated  from  a  uniform 
distribution ,  they  should  inherently  have  a  white  PSD. 

Also,  since  the  signal  being  passed  through  the  linear 
filter  is  digital,  then  this  linear  filter  must  also  be 
digital.  This  means  that  its  transfer  function  must  be 


represented  in  the  Z  domain.  To  do  this,  the  transfer 
function  must  first  be  characterized  in  the  analog  domain 
and  then  transformed  into  the  Z  domain.  To  describe  the 


analog  filter  transfer  function  the  proper  filter 
coefficients  must  be  found.  Since  the  desired  PSD  of  the 
intensity  fluctuations,  S  (w),  resembles  a  Chebychev 
filter,  then,  from  Eq .  (1)  the  filter  transfer  function  must 
also  resemble  a  Chebychev  filter.  To  properly  characterize 
this  filter,  its  order,  cut-off  and  passband  attenuation 
magnitudes  must  be  specified. 

The  order  of  the  filter  determines  how  sharp  the  roll 
off  rate  will  be  in  its  transition  band.  The  higher  the 
order,  the  sharper  the  roll  off.  By  knowing  the  various 
cutoff  frequency  and  attenuation  specifications,  the  filter 
order  can  be  computed.  To  determine  the  order,  the 
following  relation  (10:90)  can  be  used 


n  =  cosh"1  ((10’1Arain-  1 )/( 10‘  Utnax-  i))^ 

cosh*'  (  <*)  s/  ^c ) 


where 


n  =  order  of  the  filter 

A  .  =  attenuation  (db)  of  frequencies  in  the  stopband 

Amin  =  magnitude  (dB)  of  ripple  in  the  passband 
"«s  =  fre9uency  at  Amjn 

=  maximum  frequency  which  will  pass  through  the 
filter  with  a  minimum  of  attenuation 


and  A  .  „  are  constants,  while 
max  min 

which  are  functions  of  the  altitude 


w  and  w_  are  variables 

S  C 

,  propagation  path 


Chebychev  filter  orders  for  a  variety  of  optical  propagation 

environments . 


Aircraft 
Speed  (knts) 

Altitude 

(ft) 

Propagation 

Path  Length  (m) 

Filter  Order 

300 

25,000 

80, 000 

2.9639 

350 

30,000 

50,000 

2.9861 

400 

30, 000 

50, 000 

2.9473 

450 

30,000 

75,000 

2.9178 

length,  and  aircraft  speed.  Table  II  presents  a  few  values 
of  n  for  a  variety  of  environmental  parameters.  As  can  be 
seen,  the  next  highest  integer  value  under  the  heading 
"Filter  Order"  is  3  .  Therefore,  the  order  of  this  linear 
filter  can  be  approximated  to  be  3  ♦ 

A  normalized  3rd  order  analog  Chebychev  filter  will 
always  have  a  transfer  function  as  follows 

1 

H(s)  =  — r - - -  (2) 

s'3  +  1. 62933*  +  2.0773s  +  1.1516 

To  convert  this  normalized  filter  transfer  function  into  any 
desired  filter  transfer  function,  the  variable  s,  in 
Eq .  (2),  must  be  divided  by  the  frequency  at  the  desired 
filter's  3dB  point.  The  result  is 


(s/w)3  +  I.6293(s/W)2  +  2.0773(S/W)  +  1.1516 

where 

W  =  3  dB  bandwidth 


Thus,  Eq .  (3)  is  the  transfer  function  of  the  desired  analog 
filter . 

The  transformation  from  the  analog  to  the  digital 
domain  is  not  a  difficult  one.  This  conversion  can  be  done 
by  using  the  bilinear  transformation  technique.  To  obtain 
H(z),  the  following  substitution  must  be  made  into  Eq .  (3) 


Let 

where 

t  =  sampling  time 


(4) 


The  result  of  this  substitution  is 


H  ( z )  = - - - ±-3 - (5) 

/  2( z- 1 )  V  +  1.6293/  2(2-1)  r  +  2.077 U  2(z-1)  \+  1.1516 
Vt  w  (z+1 )/  \tu)  (z+1 )/  \tw  (z+1V 


Substituting  Eq.  (4)  into  Eq .  (3),  however,  does  result  in  a 
distortion  of  the  frequency  axis.  To  compensate  for  this 
distortion,  the  3  dB  frequency  term  present  in  Eq .  (5)  must 


be  warped.  Using  the  following  substitution,  this 
distortion  can  be  minimized.  Let 


W  =  2  Tan”  (  w0 1/2) 

l 


where 


t  =  sampling  time 

u>0  =  frequency  at  3  dB  point  of  the  digital  filter 


Where 


Z(z)  =  H ( z )  X ( z ) 


Thus 

H  (  z )  =  Z  (  z )  = _ _ _ _ 

TTz T  ( 1+b+C+D)z3+(C-3-B  +  3D)z2  +  (5-B-C-3D)z  +  (B-l-C+D) 


Multiplying  by 


and  rewriting 


Z (  z )  [(  1+B+C+D)  +  (C-3-B+3L-)z_1  +  (3-B-C-3D)z"2  +  (B-1-C+D)z‘3] 

=  X( z )  |a  +  3Az”^ ♦3Az“2+Az-3 ] 


(10) 


Taking  the  inverse  Z  transform,  Eq  .  (10)  may  be  written  as: 

z(n)  =  1  fkx(n)  +  3Ax(n-1)  +  3Ax(n-2)  +  Ax(n-3) 

UB+C+Z  1 

-  (C-3-B+3D)z(n-1)  -  (3-B-C-3D)z(n-2) 

-  (b-1  -C+B)z(n-3)]  (ID 

Equation  (11)  is  written  in  a  form  which  lends  itself 
extremely  well  to  computer  simulation.  As  can  be  seen,  z(n) 
is  a  function  of  past  and  present  inputs  as  well  as  past 
outputs.  Since  z(n)  is  a  series  of  random  variables,  they 
have  a  set  of  statistics  associated  with  them.  Also,  since 
the  random  variables  in  x(n)  are  zero  mean  gaussian  random 
variables,  and  the  transformation  is  linear;  then  z(n)  must 
also  be  a  sequence  of  zero  mean  gaussian  random  variables. 
However,  its  variance  will  not  be  the  same  as  that  of  x(n). 
Instead  it  will  be 

a 

al=  _!_  S  ty(n)]2  (12) 

y  a- i  n=0 


In  addition,  the  PSD  of  z(n)  will  not  be  white.  Since  the 
filter  transfer  function  is  equal  to  the  square  root  of 
S  («),  then  from  Eq .  (1),  z(n)  should  have  the  power 
spectrum  of  the  intensity  fluctuations. 

The  nonlinear  device  shown  in  Figure  13  must  now  be 
specified.  According  to  Gujar  and  Kavanaugh,  this  nonlinear 
device  is  really  a  one-to-one  transformation  between  z  and 
y.  Consequently,  for  any  point  (z,y)  in  the  transformation, 
the  probability  of  z  being  in  the  range  z-Sz  to  z  is  equal 
to  the  probability  that  y  is  in  the  corresponding  range  of 
y-jy  to  y.  That  is 

Pz(z-az/2)  1 6  Z  |  =  p  (y-ay/2)  layl  (13) 

Since  the  probability  distributions  of  z  and  y  are  known, 
then  the  only  variables  which  need  to  be  evaluated  are  bz 
and  4y.  To  solve  for  these  parameters,  a  small  value  is 
assigned  to  bz.  For  example,  let  bz-  .0001  .  To  determine 
the  corresponding  ay,  an  initial  z  and  y  value  must  be 
chosen  such  that 


Z  <  z  )  -  Py(  Y  <  y  )  ( H) 

Initial  values  of  z  and  y  are  chosen  such  that  the  upper 

limits  (i.e.  P  (  Z  <  z  )  =  P  (  Y  <  y  )  =  .999  )  of  Eq .  ( 1 4) 
z  y 

are  met.  Once  these  values  for  z  and  y,  and  Jx=  .0001  are 
substituted  into  Eq .  (11),  ay  can  be  determined.  Once  by  is 
calculated,  then  for  any  value  of  z,  the  corresponding  value 
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of  y  can  be  found.  Thus,  the  output  sequence,  y(n),  will  have  a 
lognormal  distribution. 

This  y(n)  sequence  is  then  designated  as  t(n)  and  is  shown 
in  Figure  10.  As  shown,  it  will  be  multiplied  by  the  binary 
signal,  s(n),  and  form  a  new  signal  s^(n)  which  will  have  an 
intensity  which  fluctuates  lognormally.  It  should  be  noted  that 
there  cannot  be  a  one-to-one  mul tipi ication  between  t(n)  and 
s(n).  S(n)  is  being  transmitted  at  a  rate  of  20,000  bits  per 
second  (or  1  bit  every  .05  msec),  while  the  "transmission"  of  the 
turbulence  factors,  t(n),  depends  on  the  sampling  time. 

Typically,  this  sampling  time  is  greater  than  .05  msec.  This 
means  that  the  same  turbulent  factor  number,  t(n),  might  be 
multiplied  by  more  than  one  signal  bit.  For  the  purposes  of  this 
simulation,  the  number  of  information  bits  that  will  be 
multiplied,  must  be  an  integer  value  (bit  splitting  will  not  be 
permitted).  For  this  to  occur,  the  digital  sampling  time  must  be 
such  that  when  multiplied  by  the  bit  rate,  an  integer  product 
will  result.  For  example,  assume  that  the  digital  sampling  time 
is:  t=.14  msec,  then 

20,000  bps  x  .14msec  =  2.8  (15) 

Since  2.8  is  not  an  integer,  then  the  sampling  time  must  be 
increased  so  that  the  next  highest  integer,  3>  is  obtained  in 
Eq.  (15).  To  obtain  this  value,  the  new  sampling  time  must  be 

t  =  3/20,000  bps 


.15  msec 


Thus,  for  every  turbulent  factor  produced,  it  will  be  multiplied 
by  3  consecutive  signal  bits. 


To  determine  if  the  Gujar  Kavanaugh  method  does  indeed 
produce  the  desired  PSD  at  the  output  of  the  non-linear  device, 
the  Fast  Fourier  Transform  (FFT)  can  be  used.  By  taking  the  FFT 
of  each  sequence,  z(n)  and  y(n),  and  then  comparing  them;  a 
determination  can  be  made  as  to  their  similarity.  Figure  14 
illustrates  the  FFT  for  a  sample  sequence.  Note  that  the  PSD  of 
the  sample  digital  sequence  repeats  itself  every  2000  Hz.  This 
is  because  the  sampling  time  is  .5  msec,  or  1/2000  Hz.  From  this 
figure,  it  can  be  seen  that  one  PSD  is  indeed  a  scaled  version  of 
the  other,  and  both  have  identical  frequency  components. 
Therefore,  it  can  be  concluded  that  the  Gujar  Kavanaugh  method 
will  work  for  this  simulation. 

Background  Radiation 

For  the  purposes  of  this  thesis,  it  will  be  assumed  that  the 
major  source  of  background  radiation  will  be  due  to  the  sunlit 
sky.  Typically,  the  strength  of  the  background  radiation  is 
described  in  terms  of  its  spectral  radiance,  N(A).  The  spectral 
radiance  can  be  defined  as:  the  radiant  power  in  a  solid  angle 
per  unit  projected  area  of  a  source  on  a  hemisphere  (11:113)* 
Accurate  values  for  N(  A  )  cannot 
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Figure  14.  FFT  of  the  sequences  z(n)  and 
y(n) . 


usually  be  computed  from  a  general  mathematical  expression. 
Instead,  they  must  be  measured  experimentally.  This  should 
make  sense  since  N(  X  )  is  a  function  of  many  environmental 
parameters  such  as  altitude,  cloud  cover,  temperature,  and 
source-receiver  geometry.  For  the  purpose  of  this  thesis, 
N(  X)  will  be  chosen  such  that  it  represents  somewhat  of  a 
worst  case  condition.  Once  N(  X  )  has  been  chosen,  then  the 
average  background  radiation  power  can  be  computed  (11:118) 
from 


Pb  =  Tr  N(  X  )  n  B  Ar 

where 

T  =  receiver  optics  transmissivity 
q  =  field  of  view  (steradians) 

B  =  Optical  bandwidth  (  &  )  2 

A  =  area  of  receiver  aperture  (rn  ) 

For  the  optical  communication  link  used  in  this  thesis: 

N (  X  )  =  .015  W/(m2  sr  l) 


Therefore 


fi  =  3- 8^6x10"_5,  sr 
Ar  =  1.2667x10  '  nr 
B  s  100  A 

Pb  =  (.8)(.015H3.848x10_5)(  100)(  1.2667x10"2) 
=  5.85x10”^  Watts 

This  value  represents  the  total  amount  of  power  impinging  on 
the  detector  in  the  absence  of  any  information  signal. 


Optical  Beam  Divergence 

As  the  beam  transmission  path  increases,  the  optical 
beam  begins  to  diverge.  Divergence  can  be  defined  as  the 
spreading  of  an  optical  beam  as  its  propagation  path 
increases.  As  a  result,  the  amount  of  power,  in  the 
receiver  plane,  is  smeared  over  an  area  much  larger  than 
that  of  the  transmitter  aperture.  From  Figure  15,  it  can  be 
seen  that  the  amount  of  power  received  (Pr),  is  proportional 
to  the  power  collected  over  the  total  power  present  in  the 
receiver  plane. 

This  relationship  can  be  written  as 

Pr  *  Pt  Ta  Tt  Tr  ■  <>9> 

where  «2  l2 


P.  =  transmitter  power 
A^x  =  area  of  receiver 
Lr  =  propagation  path  length 
x  =  atmospheric  transmissivity 
r  =  receiver  transmissivity 
rl  =  transmitter  transmissivity 
=  field  of  view 
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Receiver  plane 
Area  =  rr(6R )  ^ 

A 

Receiver  antenna 
Area  =  rr 

A 


Figure  15.  Divergence  of  the  optical  beam  as  the 
propagation  path  length  increases. 

The  transmissivity  factors  present  in  Eq .  (19)  represent  the 
efficiency  of  transmission  through  the  desired  medium.  The 
atmospheric  transmissivity  factor  (T_)  determines  the  amount 
of  power  which  is  able  to  pass  through  the  atmosphere 
unattenuated.  This  factor  can  be  written  (12:61)  as  follows 

T  =  exp(-  a  0L )  (20) 

3  a 

where 

aa =  extinction  coefficient 
0  =  aerosol  scaling  factor 
L  =  propagation  path  length  (km) 

From  LOWTRAN  5  computer  simulation  plots  (12:24,31)*  the 
extinction  coefficient  and  aerosol  scaling  factors  can  be 
written  as  follows: 

For  X  =  .904  m  a.  =  .45  , 

3 

0  <  altitude  <  4  km: 

8  =  .2-. 19/altitude\  (21 ) 

V  TOMT'V 
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4  km  <  altitude  <  7  km: 


B  =  .01-. 004  /altitude-4000\ 
\  3000  / 


7  km  <  altitude  <  16  km: 

5=  .006-. 0056  /al ti tude-7000\ 

\  90^  / 

+  . 0056/altitude-7000\  /altitude-7000  -1\ 
\  9000  /  \  9000  / 


As  can  be  seen  above,  the  aerosol  scaling  coefficient  is  a 
function  of  the  altitude.  The  transmitter  and  receiver 
transmissivity  factors  describe  the  percentage  of  the 
optical  signal  power  which  is  passed  through  both  the 
transmitter  and  receiver,  respectively. 


Simulation  of  the  Receiver 

Once  the  corrupted  signal's  field  is  incident  upon  the 
photodetector,  an  electrical  current  is  produced.  While 
passing  through  the  receiver,  this  current  undergoes 
additional  corruptive  effects  (i.e.  dark  current,  shot  and 
thermal  noise).  This  signal  is  then  fed  into  a  decision 
making  device  called  a  threshold  detector.  The  threshold 
detector  compares  the  level  of  the  signal  with  a 
predetermined  threshold  value.  If  the  signal  is  equal  to  or 
greater  than  this  threshold  value,  then  the  receiver  decides 
that  a  1  was  transmitted;  otherwise,  a  decision  is  made  that 
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a  0  was  sent. 


Before  developing  the  mathematical  model  for  this  type 
of  receiver,  the  following  photodetector  design  criteria 
will  be  assumed: 

1)  The  detector  will  be  thermal  noise  limited. 

2)  A  DC  current  block  will  be  used. 

3)  The  load  resistance  (R^)  will  be  equal  to 
30  ohms. 

4)  The  receiver  has  an  internal  gain  of  40. 

Since  the  dark  current  is  strictly  DC,  the  DC  current  block 
filters  it  out  from  the  average  current.  Therefore,  for  the 
purpose  of  this  simulation,  its  effects  can  be  neglected. 

but,  there  are  two  other  types  of  noise  which  must  also 
be  dealt  with.  They  are,  the  shot  noise  and  the  thermal 
noise.  If  we  assume  that  the  photodetector  is  thermal  noise 
limited,  then  we  can  overlook  the  shot  noise  due  to  the 
photodetector.  However,  we  cannot  ignore  the  shot  noise  due 
to  the  constant  background  radiation. 

Typically,  both  the  thermal  and  background  noise 
manifest  themselves  as  unwanted  current  components.  These 
extraneous  currents  can  be  modeled  as  random  variables  which 
get  added  directly  to  the  received  signal  current. 
Characteristically,  the  thermal  noise  current  can  be  modeled 
as  a  zero  mean  gaussian  random  variable,  while  the 
background  noise  current  has  been  found  to  exhibit  Poisson 
characteristics.  If  we  assume  that  the  background  noise 
level  is  quite  high,  such  as  daylight  conditions,  then  it 


can  be  assumed  that  the  background  noise  can  also  be  modeled 
as  a  gaussian  random  variable  (2: 222-22M) . 

Since  both  processes  have  a  zero  mean,  then  their 
variances  are  directly  proportional  to  there  respective 
PSD's.  In  fact,  their  variances  are  equal  to  their 
respective  power  levels.  Also,  because  of  the  fact  that  we 
are  using  avalanche  photodetectors  (APD's),  the  variances 
may  be  written  as  follows: 

For  thermal  noise, 

<y2  =  (2)(B)(2kT) 


where 

b 

=  electrical  bandwidth 

k 

=  Boltzmann’s  constant 

T 

=  temperature  °K 

R1 

=  load  resistance 

For  the  background  radiation, 

ffl=  (2)(B)(T)  Pb  q2  G2+'3) 
h  f 

where 

b  =  electrical  bandwidth 
V  =  efficiency 

P.  =  power  of  background  radiation 
q  =  charge  of  electron 
G  =  gain 

h  =  Planck’s  constant 
f  =  frequency  of  light 


Using  the  parameters  of  the  communication  system,  the 
respective  variances  are 


<6.62x10"34)(3.31x1014) 

=  2.804x10“15 

Now,  by  using  the  IMSL  routines,  these  normal  (0, 

2.76x10-1^)  and  normal  (0,2.804x10“^)  random  variables  can 
be  generated  and  then  added  to  the  signal  current.  The 
resulting  signal  is  then  passed  through  the  threshold 
device . 

The  threshold  value  for  this  device  must  be  chosen  such 
that  the  probability  of  making  an  incorrect  decision  in 
minimized . 

Determination  of  the  Threshold  Level 

In  order  to  determine  the  effects  of  atmospheric 
turbulence  on  an  optical  signal,  a  "baseline"  system  must 
first  be  developed.  The  design  of  this  baseline  system  will 
consist  of  a  communication  link  which  does  not  take  into 
account  the  effects  of  turbulence.  By  comparing  a 
turbulence  corrupted  signal  with  a  signal  from  this  baseline 
system,  a  comparison  can  be  drawn.  This  comparison  will 
enable  us  to  determine  how  the  bit  error  rate  is  affected  by 


turbulence . 


Since  we  are  simulating  an  on/off  keying  communication 


system,  then  the  only  two  possible  hypothesis  m^  and  m2 
corresponding  to  the  signals  entering  the  decision  device 
are 


1  : 

s  =  B  +  T 

(24) 

2  1 

s  =  n  +  B  +  T 

(25) 

where 

T  =  random  thermal  noise  current 
b  =  random  background  radiation  current 
M  =  average  signal  current 

Since  both  T  and  B  are  zero  mean  gaussian  random  variables 

Let  X  =  B  +  T 

Where  E { X }  =  0 

e{x2}=  *1  +  al 

Therefore  =  2. 804x 1 0“ 15  +  2.76x1(T15 

=  5.564x10-15 

This  allows  Eqs .  (24)  and  (25)  to  be  rewritten  as 

m 1 :  s  =  X  (26) 

m2:  s  =  m  +  X  (27) 


Now  p ( s [ m i )  =  px ( X) 

p(s |m2)  =  px (X-  M  ) 

Since  the  probability  distribution  of  X  is  gaussian,  we  get 


PXU) 


727*; 


exp/-1 

l  ToL 


(28) 


(29) 


Pv  (x  -  m  )  =  1  exp  /-I  (x-m)2\ 

lr^x  / 

The  system  probability  of  error,  P  ,  can  be  written  as 

Pe  =  P(deciding  a  1  J  0  sent)  P  ( 0 )  + 

P(deciding  a  0  |  1  sent)  P(l) 

Let  P(0)  =  P(1)  =  h 

Then 

Pe  =  ^  (  p(decidin8  a  1  |  0  sent)  + 

P(deciding  a  0  |  1  sent)) 


where 


k  =  threshold  value 


The  optimum  threshold  (11:169)  value  for  this  type  of  system 


would  be 


where 


H  -  average  current  due  to  the  signal  , 
and  the  background  radiation 


(3D 


( 


Before  a  threshold  value  can  be  assigned,  an  optimum  value 
of  the  system  probability  of  error  must  be  designated.  Due 


to  the  constraints  of  computer  time,  Pg  for  this  system  must 
be  at  least  lO-'*. 

Substituting  Eq .  (31)  into  Eq .  (30)  and  rewriting  gives 


(32) 


If,  for  example 

Q ( x )  =  ID'5 

Then  x  =  *4. 265 


Therefore 


4.265  =  U 
2*x 

Or  H  -  ( 4. 265) (2) (  5.564x10'15) 

H  =  6. 36x 1 0"7  Amps 

To  determine  the  proper  threshold  level,  we  would  use  the 
value  calculated  for  M  ,  and  then  substitute  it  into 


Eq.  (3D. 


Thus 


k=  (6 . 36x 1 0~7  ) 
2 


=  3 . 1 8 x 10~7  Amps 


Therefore,  the  threshold  level  must  be  set  at 
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3.18x10  Amps.  By  changing  the  average  current  value,  n, 
in  .] .  (32),  and  using  Eqs.  (  19),  (20),  (21  ),  (22),  and  (23) 
the  bit  error  rate  can  be  determined  for  any  altitude  and 
path  length.  Using  the  fixed  threshold  value  above, 
theoretical  bit  error  rates  for  an  increasing  propagation 
path  length,  at  any  altitude,  can  be  determined. 

Summary 

In  modeling  the  effects  of  the  turbulence,  the  Gujar 
kavanaugh  technique  was  used  to  produce  a  sequence  of 
numbers  which  have  a  lognormal  distribution  and  a  low  pass 
filter-like  PSD.  By  taking  the  FFT  of  this  turbulence 
sequence,  it  was  found  that  its  PSD  actually  resembled  that 
of  the  intensity  fluctuations.  Once  these  sequence  of 
numbers  were  generated,  they  were  multiplied  by  the  binary 
information  signal.  The  result  was  a  signal  which  faded 
with  time.  After  undergoing  the  multiplicative  effects  of 
turbulence,  the  effects  on  the  signal  due  the  background 
radiation,  were  evaluated.  It  was  found  that  the  power 
associated  with  the  background  radiation  was  constant 
and  a  function  of  its  spectral  radiance.  The  receiver  was 
assumed  to  be  thermal  noise  limited  and  the  shot  noise  due 
to  the  background  radiation  was  assumed  to  be  normally 
distributed.  Thus,  gaussian  statistics  could  be  used  in 
determining  the  bit  error  rate  and  the  optimum  threshold 


IV.  Results  of  the  Simulation 


Introduction 

This  chapter  examines  the  performance  of  the  on/off 
keying  air-to-air  optical  communication  link.  The 
performance  results  for  this  link  have  been  presented  in 
terms  of  the  bit  error  rate  plotted  against  the  propagation 
path  lengths. 

In  determining  the  performance  of  this  communication 
system,  a  total  of  nine  different  links  were  chosen  to  be 
evaluated.  The  performance  parameters  associated  with  these 
links  are  listed  in  Table  III.  On  the  whole,  these  links 
were  chosen  such  that  they  represented  potential  operating 
environments  of  the  larger  operational  aircraft  in  the  Air 
Force  inventory.  By  varying  the  altitude,  aircraft  speed, 
and  data  rate  a  basic  understanding  of  how  these  parameters 
can  affect  the  performance  of  the  link  can  be  achieved. 

The  next  logical  step  in  the  performance  determination 
is  to  develop  a  baseline  linK.  For  the  purpose  of  this 
thesis,  the  baseline  link  will  be  assumed  to  operate  under 
the  identical  conditions  of  the  simulated  air-to-air  link; 
except  that  it  will  not  be  influenced  by  the  effects  of 
turbulence.  By  comparing  the  air-to-air  link  with  the 
baseline  link,  an  understanding  of  how  atmospheric 
turbulence  affects  these  optical  links  can  be  achieved. 
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Table  III 


Parameters  of  the  communication  links  used 


Altitude 

(ft) 

True  Aircraft  Speed 
(knots ) 

Data  Rate 
(bps) 

10, 000 

425 

20, 000 

10,000 

425 

40, 000 

10, 000 

475 

20,000 

20,000 

425 

20, 000 

20, 000 

425 

40, 000 

20,000 

475 

20,000 

30,000 

425 

20, 000 

30, 000 

425 

40, 000 

30, 000 

475 

20, 000 

Table  IV 

Transmitter  and  Receiver  Characteristics 


Transmitter  Output  Power:  50  W 
Transmitter  Optics  Efficiency:  0.9 

Transmitter  and  Receiver  Aperture  Diameters:  0.127  m 
Receiver  Optics  Efficiency:  0.8  ^ 

Receiver  Field  of  View:  3*848x10~'>  sr 
Receiver  Gain :  MO 


The  equipment  parameters  used  in  this  simulation  are 
listed  in  Table  IV.  If  desired,  the  choice  of  these 
parameters  could  be  purely  arbitrary.  Changing  any  of  the 
parameter  values  doesn't  cause  any  problems  with  the 
simulation  software. 
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Originally,  the  software  for  this  simulation  was 
written  on  AFIT’s  VAX  computer.  However,  the  sheer  number 
of  calculations  and  IMSL  subroutine  calls  proved  to  be  much 
more  than  the  VAX  could  efficiently  handle.  As  a  result, 
Aeronautical  Systems  Division's  (ASD)  CYBER  computer  system 
had  to  be  used.  This  system  proved  to  be  more  capable  of 
running  the  simulation. 

Performance  of  the  Baseline  System 

By  using  the  bit  error  rate  analysis  described  in 
Chapter  III,  baseline  link  performance  curves  were 
generated.  Although  these  curves  do  not  include  the  effects 
of  turbulence;  thermal  current  and  shot  noise  due  to  the 
background  radiation  are  accounted  for. 

The  baseline  link  performance  curves  are  plotted  in 
Figures  16,  17,  and  18.  As  shown,  the  bit  error  rate  rises, 

_5 

from  5x10  to  . 5  ,  as  the  propagation  path  length 
increases.  This  "fluid"  bit  error  rate  is  due  to  the  fact 
that  as  the  path  length  increases,  the  received  signal 
strength  decreases.  As  a  result,  the  thermal  current  and 
the  shot  noise,  due  to  the  background  radiation,  begin  to 
dominate  over  the  signal. 

Since  our  threshold  levels  have  been  designed  for  a 

-5 

minimum  bit  error  rate  of  10  ,  it  can  be  said  that  the 

system  has  a  certain  "reliability  range."  This  reliability 
range  corresponds  to  the  maximum  range  that  the  link  can 
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operate  while  still  maintaining  a  bit  error  rate  of  at  least 
10  .  From  Figures  16,  17,  and  18,  it  can  be  seen  that 

these  maximum  reliable  ranges  are  95  km,  208  km,  and  239  km, 
respectively . 

Although  the  received  power  of  the  signal  is  the 
dominant  factor  in  determining  the  bit  error  rate,  this 
power  is  a  function  of  many  factors.  Factors  such  as:  path 
length,  field  of  view,  and  receiver  optics  play  key  roles  in 
determining  what  the  received  power  levels  will  be. 
however,  the  absorption  and  scattering  coefficients  are  the 
factors  that  allow  each  link  to  obtain  longer  maximum 
reliable  ranges  with  an  increase  in  the  altitude.  These 
coefficients  determine  the  amount  of  transmitted  signal 
power  which  passes  through  the  atmosphere,  unattenuated,  to 
the  receiver.  The  higher  the  value  of  these  coefficients, 
the  larger  the  amount  of  received  signal  power.  As  the 
altitude  is  increased,  these  values  typically  increase. 

By  looking  at  these  performance  curves  it  becomes 
obvious  that  the  maximum  reliable  range  gets  longer  as  the 
altitude  is  increased.  This  is  primarily  due  to  the  fact 
that  there  is  less  of  a  concentration  of  aerosols  (i.e. 
dust,  smoke,  smog,  etc...)  in  the  higher  altitudes.  Thus, 
there  is  less  scattering  of  the  optical  beam. 


The  computer  simulation  was  run,  and  the  necessary 
performance  data  points  were  generated  for  the  nine  links 
shown  in  Table  III.  The  results  of  these  runs  were  broken 
up  into  three  altitude  groups;  namely,  10000  ft,  20000  ft, 
and  30000  ft.  In  turn,  each  of  these  groups  were  broken 
down  into  three  subgroups;  where  each  subgroup  consists  of 
links  with  varying  aircraft  speeds  and  data  rates. 

The  results  of  the  computer  simulation  for  the 
10,000  ft  altitude  group  can  be  seen  in  Figures  19,  20,  and 
21,  and  the  results  for  the  20,000  ft  links  are  shown  in 
Figures  22,  23,  and  24.  Figures  25,  26,  and  27  present 
those  results  for  30,000  ft.  The  data  points  used  to 
generate  these  performance  plots  are  listed  in  Appendix  A. 

By  comparing  Figures  19  -  27  with  the  baseline  links 
shown  in  Figures  16,  17,  and  18,  the  effects  of  turbulence 
can  clearly  be  seen.  As  shown,  atmospheric  turbulence 
severely  degrades  the  performance  of  the  optical  link. 

Table  V  compares  the  maximum  reliable  communication  ranges 
of  the  turbulence  induced  link  vs.  those  of  the  baseline 
link.  As  shown,  the  maximum  reliable  ranges  are  degraded 
anywhere  from  50  to  66  percent;  depending  upon  the  link 
analyzed.  This  amount  of  degradation  could  have  a  severe 
impact  on  the  operational  environment  of  any  optical 
communication  system. 
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Figure  21 


Turbulence  affected  optical  link 
Altitude:  10,000  ft 

Aircraft  speed:  425  knts 
Data  rate:  40,000  bps 
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22.  Turbulence  affected  optical  link 
Altitude:  20,000  ft 
Aircraft  speed:  425  knts 
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Figure  23.  Turbulence  affected  optical  link 
Altitude:  20,000  ft 

Aircraft  speed:  ^75  knts 
Data  rate:  20,000  bps 
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Figure  25-  Turbulence  affected  optical  links. 
Altitude:  30,000  ft 

Aircraft  speed:  425  knts 
Data  rate:  20,000  bps 
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Figure  21.  Turbulence  affected  optical  link. 
Altitude:  30,000  ft 

Aircraft  speed:  425  knts 
Data  rate:  40,000  bps 


Table  V 


Maximum  reliable  ranges  for  the  baseline  link 
and  the  turbulence  affected  links 


Further  analysis  of  the  performance  curves,  and  the 
data  in  Appendix  A,  of  each  altitude  groMD,  yielded  some 
interesting  results.  Namely,  the  plots  of  each  altitude 
group  show  that  as  the  aircraft  speed  is  increased,  from 
425  knots  to  475  knots,  the  bit  error  rates  didn't  appear  to 
change  significantly.  This  could  be  due  to  the  fact  that 
the  change  in  velocity  (50  knots)  was  relatively  small. 

However,  the  fact  remains  that  the  speed  of  the 
aircraft  doesn’t  affect  the  strength  of  the  turbulence. 
Instead,  variations  in  speed  will  only  effect  the  frequency 
content  of  the  intensity  fluctuation  power  spectral  density 
(IFPSD).  From  Eq .  (20)  in  Chapter  II,  it  was  shown  that  the 
highest  frequency  component,  f  ,  present  i>~,  the  IFPSD  can  be 


expressed  as 


where 


ws  =  aircraft  speed 
L  =  propagation  path  length 

From  this  relation,  it  can  be  seen  that  for  equivalent  path 
lengths,  the  IFPSD  contains  higher  frequency  components  for 
faster  aircraft  speeds.  This  means  that  the  optical  beam's 
intensity  fluctuations  will  vary  more  rapidly  as  the  speed 
of  the  aircraft  is  increased. 

An  example  of  this  is  shown  in  Figure  28.  As  shown, 
the  first  32  msecs  of  the  200  knot  turbulence  factor 
sequences  have  been  compressed  (and  are  equivalent)  to  the 
first  11  msecs  of  the  600  knot  sequences.  This  clearly 
shows  that  as  the  speed  is  increased,  the  turbulence  factors 
vary  more  quickly.  However,  their  maximum  and  minimum 
values  are  not  affected  by  the  aircraft's  speed.  From  this 
figure,  it  can  be  determined  that  the  variance  of  the  log 
amplitude  is  responsible  for  inducing  larger  fluctuations  of 
turbulence  factor  values. 

A  sample  turbulence  factor  sequence  for  425  knots  and 
475  knots  is  shown  in  Figure  29.  As  shown,  the  sequence  at 
475  knots  is  slightly  more  compressed  than  the  one  at 
425  knots.  This  close  similarity  in  these  sequences  is  the 
reason  for  similar  bit  error  rates  at  different  aircraft 
speeds.  One  could  then  reasonably  assume  that  as  the 
difference  in  aircraft  speeds  got  larger  the  bit  error  rates 
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Aircraft  Speed:  200  knots 
Variance:  .0976 
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Aircraft  Speed:  425  knots 
Variance :  .1469 
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b.  Aircraft  Speed:  475  knots 
Variance :  . 1469 


Figure  29.  Representation  of  the  turbulence  factor 
sequences  for  an  altitude  of  20,000  ft 
and  a  70  km  path  length. 


could  be  expected  to  differ  by  a  more  significant  amount. 
However,  one  must  remember  that  the  larger  and  more  typical 
aircraft  in  the  Air  Force  fleet  tend  to  cruise  at  speeds 


Thus,  the  variance  of  the  log  amplitude  is  the  dominant 
factor  in  determining  the  values  of  the  turbulence  factors. 
From  the  data  in  Appendix  A,  it  can  also  be  shown  that  as 
the  variance  gets  ranging  from  375-500  knots;  a  difference 
of  only  125  knots.  Therefore,  it  can  be  concluded  that 
increasing  or  decreasing  the  speed  of  the  aircraft  will  not 
have  a  dramatic  impact  on  the  link  bit  error  rates. 

Thus,  the  variance  of  the  log  amplitude  is  the  dominant 
factor  in  determining  the  values  of  the  turbulence  factors. 
From  the  data  in  Appendix  A,  it  can  also  be  determined  that 
as  the  variance  gets  larger,  the  "strength"  of  the 
turbulence  increases,  and  the  bit  error  rate  rises.  From 
Chapter  II  it  was  shown  that  the  variance  could  be  expressed 
as : 


<C  = 


124  Cn  (k) 


7/6 


(L) 


11/6 


where 

k  ^  =  wavenumber 

Cn^  =  refractive  index  structure  constant 
L  =  propagation  path  length 


Since  "k"  is  a  constant  throughout  the  entire  simulation, 

2 

the  variance  becomes  a  function  of  both  Cn  and  the  path 
length.  Values  of  Cn2  for  altitudes  of  10000  ft,  20000  ft, 
and  30000  ft  have  been  calculated  to  be  b.MxlO'1^, 
1.6l7x10-1^,  and  1.947x10”^®,  respectively.  By  comparing 


the  differences  between  Cn  for  each  of  the  altitudes 

analyzed,  it  is  easy  to  see  that  they  differ  by  multiples  of 

an  order  of  magnitude.  Since  the  propagation  path  lengths, 

for  the  various  links,  differ  from  each  other  by  no  more 

2 

than  factors  of  2  and  3,  Cn  can  be  considered  to  be  the 

dominant  factor  in  determining  the  value  of  the  variance. 

2 

However,  the  value  of  Cn  is  determined  by  the  altitude  of 

the  aircraft.  Therefore,  it  is  the  altitude  which  plays  the 

key  role  in  determining  the  range  of  reliable  communications 

and  the  bit  error  rate.  The  higher  the  altitude,  the  longer 

the  reliable  range  becomes. 

However,  one  cannot  keep  increasing  the  altitude 

indefinitely.  For  example,  Figure  4  (page  16),  illustrates 

that  thermal  inversions  within  the  tropopause  which  occur  at 

2 

approximately  35,000  ft  cause  the  values  of  Cn  to  increase. 
As  a  result,  the  variance  increases  and  the  maximum  reliable 
ranges  decrease.  Therefore,  for  a  given  optical 
communication  system,  a  prudent  choice  of  altitudes  could 
minimize  the  bit  error  rates,  while  at  the  same  time 
maximizing  the  reliable  communication  ranges. 

Increasing  the  data  rate,  proved  to  consistently  raise 
the  bit  error  rate  in  all  of  the  links  analyzed.  Examining 
this  situation  further,  it  was  found  that  increasing  the 
data  rate  forced  more  bits  to  be  multiplied  by  a  singular 
turbulence  factor.  For  example,  if  we  had  1  turbulence 
factor  for  every  5  signal  bits,  then  doubling  the  data  rate 


would  yield  1  turbulence  factor  for  every  10  signal  bits. 

If  there  happened  to  be  a  prolonged  fade,  for  example,  3  or 
4  sequential  turbulence  factors  having  very  small  values, 
then  there  is  a  potential  to  "wipe-out”  a  string  of  40 
sequential  signal  bits. 

Summary 

In  conclusion,  the  computer  simulation  was  run,  and  the 
results  were  compared  with  those  of  the  baseline  link. 

Also,  assessments  were  made  as  to  how  atmospheric  turbulence 
affects  the  air-to-air  optical  link.  The  next  chapter, 
Chapter  V,  summarizes  the  conclusions  reached  from  this 
project,  and  presents  recommendations  for  future  work. 


V.  Conclusions  and  Recommendations 


Conclusions 

The  goal  of  this  project  was  to  determine  the 
performance  of  an  air-to-air  optical  communication  link,  and 
to  evaluate  the  effects  of  atmospheric  turbulence.  A  major 
part  of  this  task  was  the  development  of  the  computer 
simulation  which  determined  the  bit  error  rates  from  various 
optical  links.  Using  this  simulation,  performance  results 
were  obtained  for  nine  different  turbulence  influenced 
links.  By  comparing  the  data  from  these  links  with  the 
baseline  links,  it  was  possible  to  evaluate  the  effects  of 
turbulence.  The  results  of  both  the  baseline  link  and  the 
computer  simulation  links  have  been  presented  in  the  form  of 
performance  curves.  A  summary  of  the  conclusions  reached 
are  presented  as  follows 

1.  By  comparing  the  baseline  link  with  the  turbulence 
influenced  links  it  was  found  that  atmospheric  turbulence 
severely  degraded  the  performance  of  the  analyzed 
communication  link.  Losses  of  from  50  to  66  percent,  over 
the  theoretical  maximum  in  reliable  communication  ranges 
were  not  uncommon. 

2.  Although  the  speed  of  the  communicating  aircraft  was 
only  increased  by  50  knots,  the  differences  in  bit  error 
rates  (at  comparable  altitudes)  was  insignificant.  This  was 


due  to  the  fact  that  the  aircraft's  speed  only  affects  the 
rapidness  of  the  intensity  fluctuations  and  not  the 
strength.  Larger  differences  in  aircraft  speed  might 
produce  more  significant  results. 

3-  The  variance  of  the  log  amplitude  was  found  to  be  the 
dominant  factor  in  determining  the  strength  of  the 
turbulence.  Larger  values  of  the  variance  tended  to 

2 

increase  the  link  bit  error  rate.  It  was  shown  that  Cn 

plays  a  key  role  in  the  determination  of  the  value  of  the 

2 

variance.  Since  Cn  is  directly  related  to  the  altitude,  a 
prudent  choice  of  altitudes  could  provide  longer  reliable 
communication  ranges. 

M.  Increasing  the  data  rates  of  the  links  tended  to 
increase  the  bit  error  rates.  This  was  due  to  more  signal 
bits  being  multiplied  by  fewer  turbulence  factors.  If  the 
data  rate  were  slowed  down,  then  (for  the  purposes  of  this 
simulation)  there  would  be  fewer  signal  bits  being  affected 
by  a  singular  turbulence  factor.  Ideally,  one  would  want 
the  data  rate  at  such  a  speed  that  multiple  turbulence 
factors  would  affect  a  single  signal  bit.  Then,  using  an 
integrate  and  dump  type  receiver,  the  signal  bit  could  be 
recovered.  Slowing  the  data  rate,  in  this  way,  allows  us  to 
average  the  effects  of  the  turbulence  factors.  However, 
there  is  a  limit  as  to  how  slow  a  potential  user  would  like 
his  data  rate  to  be. 


5.  It  was  also  found  that  the  saturation  region  was  never 
reached  until  bit  error  rates  of  approximately  .3  -  .35  . 
This  means  that  within  reliable  communication  ranges 
lognormal  statistics  can  always  be  used  to  evaluate  the 
effects  of  turbulence  for  a  communication  link. 

Recommendations 

Several  follow-on  efforts  have  become  evident  as  a 
result  of  the  work  performed  during  the  period  of  this 
thesis.  The  following  recommendations  are  submitted 

1.  The  analysis  of  the  communication  link  was  based  on  an 
on/off  keying  system.  Since  these  types  of  systems  aren't 
in  common  use  today,  the  simulation  could  be  modified  to 
include  various  other  types  of  modulation  schemes.  Once 
these  other  schemes  have  been  programmed,  a  performance 
comparison  of  one  scheme  verses  another  could  be  done.  As 
result  of  these  comparisons,  a  recommendation  could  be  made 
as  to  which  modulation  scheme  provides  the  best  performance 

2.  The  present  simulation  does  not  include  any  coding  of 
the  information  signal.  An  investigation  into  the  various 
coding  and  interleaving  schemes  could  be  performed.  This 
might  give  some  insight  into  which  coding  schemes  might 
prove  to  be  the  most  effective  against  turbulence. 


3.  This  simulation  also  assumes  a  simple  threshold  detector 
as  the  decision  making  device.  It  might  prove  to  be  more 
realistic  if  a  matched  filter  device  were  used  instead. 

This  could  be  done  by  modifying  the  computer  simulation. 
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FORTRAN  Program  Commputer  Code  Listing 


PROGRAM  SINT 

*  PROGRAM  SINT  WRITTEN  BY  *  JAY  N.  KANAVOS,  CAF'T,  USAF  * 

*  THIS  PROGRAM  WILL  COMPUTE  THE  BIT  ERROR  RATE  FOR  AN  * 

*  OPTICAL.  AIR-TO-AIR  LINN  GIVEN  THE  AIRSPEED*  ALTITUDE,  * 

*  PROPAGATION  PATH  LENGTH,  AND  DATA  RATE.  * 

IK###*##*#**#*#*######*##*###***##**#****#*#***###***##*#**### 

C234567890 

DOUBLE  PRECISION  SEED(2500) 

DIMENSION  D1(2),Z(21>,X<21),F(20),B(20>,Y(21> , SUM (834) 

REAL  L.AMDA 

COMMON  FILCOA, F ILCOB,FILCOC,FILCOD,ZETA , VAR ,L AMDA ,THERMN 
+  ,  SHOTNO ,  GOI  .D1 ,  RXPOW ,  ROWS  I G ,  THRES ,  RLOAD 
PRINT  *,' INPUT  AIRCRAFT  SPEED( KNTS > , ALT ( FT > , PATH  LENGTH ( KM ) 
READ  *,V,F’,A 
VN=V*. 51398 
PATH-F'flOOO. 

ALT-A*. 30473 
N=625 
C0UNT=0 . 

*  THIS  ROUTINE  CALCULATES  THE  SEED  NUMBERS  FOR  THE  IMSL 

*  SUBROUTINES. 

DO  10  J- 1,2500 

CALL  GGUBS( 12345678. DO, 2, D1 ) 

SEED< J)=IDINT<D1(1)*D1(2)*1.D8) 

10  CONTINUE 

*  THE  FOLLOWING  REPRESENTS  THE  EQUIPMENT  PARAMETERS 

EFF= • 85 

LAMDA=9040 .E-10 

CONST =6.62E-34* ( 3 .E8/LAMDA ) 

RL0AD=30. 

Q= 1 . 6E-19 
BAND=5.E6 
P0WN-5.85E-7 
GAIN=40. 

*  THERMN  a  THERMAL  NOISE  POWER 

*  SHOTNO  =  BACKGROUND  RADIATION  SHOT  NOISE  VARIANCE (POWER) 

*  THRES  =  THRESHOLD  LEVEL  OF  RECEIVER 

THERMN=2 . *BAND*2 . *1 ♦ 38E-23*300 . /RLOAD 
SH0TN0=2.*BAND*EFF*P0WN*Q**2*GAIN**(2+.3>/C0NST 
THRES* . 5*4 . 264890793924*2 . *SQRT <  SHOTNO+THERMN ) 

PI=3. 1415926 

*  BR  *  BIT  RATE 

BR*20000. 

X< 1 )=-  .5 

RXAREA=(PI/4. >*.127**2 

*  BETA  *  AEROSOL  SCALING  FACTORS 

IF< ALT. LE. 4000. ) BETA* . 2- ♦ 19* < ALT/4000 . ) 


IF (ALT  .GT .4000*  ♦AND.  ALT.LE.7000. >  BETA- ♦ 01- ♦ 004#  < ( ALT 
+-4000. )/3000.  ) 

IF ( ALT. GT. 7000.)  BETA=6 . E-3-5 . 6E-3# ( (ALT-7000. )/9000. ) 
++. 5* ( (ALT-7000. )/9000. )*(  ( (ALT-7000. )/9000. >-l.)*1.12E- 

*  ATRAN  ■  TRANSMISSIVITY  OF  ATMOSPHERE 

*  FOV  =  FIELD  OF  VIEW  (SR) 

*  RXPOW  =  RECEIVED  SIGNAL  POWER 

*  POWSIG  =  AMOUNT  OF  SIGNAL  CURRENT  IN  RECEIVER 

ATRAN=FXP ( - . 42*BET A#P ) *EFF# . 67 

FOV= ( PI/4 . )*7.E-3**2 

RXPOW- ATRAN#50 . *RXAREA/ ( F  0V*PATH**2 ) 

P0WSIG=(EFF*RXP0W*Q*GAIN/C0NST)**2 

*  THIS  SECTION  WILL  CONVERT  THE  NORMMALIZED  POWER  SPECTRUM 

*  INTO  A  DENORMALIZED  POWER  SPECTRUM 


DO  13 

J=2 , 20 

X(  J)  = 

X( J-l )+.0625 

CONTINUE 

B  ( 1  )  = 

.  4320 

B  ( 2  )  = 

.4333 

B  ( 3  )  = 

.4433 

B(4)  = 

.4467 

B(5)  = 

.45 

B  (  6 )  = 

.448 

B(7)  = 

.4467 

B(  8  )  = 

.  44 

B(9)  = 

.432 

B(  10) 

=.4133 

B  ( 1 1 ) 

=.3833 

B(  12) 

= . 3333 

B  ( 13 ) 

=  .2667 

B(  14) 

=.1833 

B(  15) 

=.0733 

B(  16 ) 

=  .065 

B(  17) 

=.0407 

B(  18) 

=  .03 

B(  19) 

=  .02 

B(  20) 

=  .015 

************************************************************ 


*  THE  VALUES  FOR  EACH  POINT  ON  THE  PSD  PLOT  ARE  * 

*  DETERMINED  HERE.  * 

*  VN  ■  WIND  SPEED  * 

*  L  =  PATH  LENGTH  * 

*  FM  *  MAXIMUM  FREQUENCY  * 

*  FO  =  A  NORMALIZATION  FREQUENCY  * 


*********************] *************************************** 

*  CN2  *  REFRACTIVE  INDEX  STRUCTURE  CONSTANT 

*  VAR  =  VARIANCE  OF  THE  LOG  AMPLITUDE 

I F < ALT . LT . 1000 . )  CN2= ( ( 5 . 3E-8-5 . 07E-8* (( ALT-10 . ) /990 . ) > 
+*4.6416)**2 


IF < ALT. GE. 1000.  .AND.  ALT .LT .8000. )  CN2=< (2.3E-9+ 

+  ( <  ALT- 1000 ♦ )/7000* )*(-1.97E-9))#4.64t6>  ##2 
IF < ALT. GE. 8000.  .AND.  ALT. LT. 10000. >  CN2=<  <3.3E-10+< (ALT- 
+8000. >/1000. >#<-3.E-ll>+<  < ALT-8000. )/1000. >*< (< ALT-8000. > 
+/1000. )-l . )#. 5*6. E-ll >#4.6416) **2 
IF ( ALT. GE. 10000.  .AND.  ALT. LE . 14000. )  CN2=( <3.3E-10+( (ALT- 
+10000. )/2000. ) #3 . 17E-9+ <  < ALT- 10000. )/2000. )♦< ( (AI  T-10000. ) 
+/2000. )-l. >*.5*<-6.34E-9) >*4.6416>*#2 
VAR- . 1 24#CN2*  <  2 . #PI/LAMDA ) **  <  7 . /6 . >  * ( PATH ) ## ( 1 1 . /6 . ) 
FM=.32#VN/(LAMDA*PATH)**< .5) 

F0=VN/SQRT(PATH*L.AHDA*2.*PI ) 

***************#***##*********#***#)^*#**#**#**#*********^ 


*  F (  )  =  FREQUENCY  VALUE  * 

*  Y(  )  =  CORRESPONDING  PSD  VALUE  * 

*  Y ( 2 1 ) =  THE  LAST  POINT  WITHIN  THE  .2  DB  OF  THE  * 

*  MAXIMUM  VALUE  OF  THE  PSD  * 


I********************************************************* 

DO  11  J=1 . 20 
F< J)=FO*EXP(X( J) > 

Y( J)=SQRT(B< J)*VAR/FO> 

11  CONTINUE 

IF  (FM.LT .F(5) )  PRINT  *, 'SORRY  SOMETHING  IS  WRONG' 
DEI..TA1=F(A)-F  <5) 

DELTA2=FM-F<5) 

DELT  A3=DELT  A2/DELTA1 
DIF1=Y(5)-Y(6) 

DZF2-DIF18DELTA3 

Y(21)=DIF2+Y(5) 

DO  12  J=1 ,20 
Z< J)=Y< J)/Y(21) 

12  CONTINUE 

*  THE  3  DB  POINT  IS  CALCULATED  HERE  * 

jQ,  j||»  ||f  j|^  j|^  ||^  j|^  ^  ^  ^  j|^  ^  ^  |||  ^  ^  ^  ^ 

DB-Y  <  21 )/2 . 

DO  33  J=l,20 
IF  (DB.GE.Y(J))  GO  TO  99 
GO  TO  33 
99  A=Y ( J) 

BZ*Y< J-l) 

C=F< J) 

D=F ( J-l ) 

GO  TO  98 
33  CONTINUE 

98  A1=BZ-A 

A2=BZ-DB 
BZ1=C-D 
BZ2=A2*BZ1/A1 
S=D+BZ2 


SAMP=  1 . / <  3 ♦ #F  <  20  > ) 

********************************************************** 

*  THIS  ROUTINE  CALCULATES  THE  NEW  SAMPLING  TIME.  * 

*  WE  NEED  A  NEW  SAMPLING  TIME  BECAUSE  THE  SCINTILLATION* 

*  MULTIPLICATIVE  NUMBERS  MUST  BE  USED  AN  INTERGER  * 

*  MULTIPLE  NUMBER  OF  TIMES.  * 

*  BR  =  BIT  RATE  <20000  B/S)  * 

GOLD=BR*SAMP 
G0LD1-AINT (GOLD) 

IF < G0LD1 . GE . 1 . )  G0LD1=G0LD1  +  1 . 

IF  <  G0LD1 . LT ♦ 1 . >  G0LD1  =  1 » 

T=G0LD1/BR 

*********************************************************** 

*  RANDOM  NUMBERS  ARE  GENERATED  AND  CHANGED  INTO  * 

*  NORMAL  DISTRIBUTED  NUMBERS.  THESE  NUMBERS  ARE  # 

*  THEN  FED  INTO  THE  FILTER,  AND  OUTPUTED.  THE  * 

*  OUTPUTED  NUMBERS  HAVE  THE  DESIRED  PSD.  * 

*********************************************************** 
********************************************************* 

*  THE  FILTER  COEFFICIENTS  ARE  DETERMINED  HERE  * 

951  S=<2./T)*TAN(T*S*PI ) 

C0EF1=. 81463413 
C0EF2=1. 41362877 
C0EF3=C0EF1 
ALPHA=COEF 1 +C0EF3 
BET  A=COEFJ *COEF3+C0EF2 
DEL.TA=C0EF2*C0EF3 
FILC0A=S**3*<T/2. >**3 
FILC0B=ALPHA*S*<T/2. ) 

FII.  C0C=BETA*(S**2>*<T/2.  )**2 
FILCOD -DELTA*  <  S**3 )*  < <T/2. >**3) 

ZETA=1 ./< 1 ,+FILCOB+FILCOC+FILCOD) 

C0UNT1=0. 

*  SUBROUTINE  'ERROR'  ACTUALLY  DETERMINES  THE  VALUES  OF 

*  OF  THE  TURBULENCE  FACTORS 

DO  1  J=1,N 

CALL  ERROR  <  ALT , J, SEED <  J ) , SEED <  J+625 ) , SEED  <  J+ 1 250 ) , 
+SEED<  J+1875) ,SUM<  J> , NR ) 

C0UNT1=C0UNT1 +NR 
COUNT =COUNT +SUM <  J ) 

1  CONTINUE 

PERROR=COUNT /COUNT 1 
SHANON= ALT/, 30473 

WR I TE  <  * , 4 ) PERROR , PATH , COUNT , SHANON , V 
4  FORMAT <  E15.8,3X,F10.2,2X,F10.2,3X,F10.1 , 2X,F5. 1 ) 

END 

*  IN  ORDFR  TO  ALLOW  FOR  EFFICIENT  USE  OF  COMPUTER  TIME, 


*  ONLY  1200  TURBULENCE  FACTORS  CAN  BE  DETERMINED  AT 

*  ONCE.  THUS,  SUBROUTINE  'ERROR'  IS  CALLED  ABOUT 
*.  750,000/1200  ■  A25  TIMES 

SUBROUTINE  ERROR < ALT , INDIA , DSEED , DSEED1 , DSEED2, 

+DSEED3 , SUM , NR  > 

DOUBLE  PRECISION  DSEED , DSEED1 , DSEED2 , DSEED3 
DIMENSION  XI <-2t 1200), Y1 (-21 1200) ,R1( 1200) , 

+R2 < 1200 ) , S IGCUR  < 1 200 ) , TOTAL  < 1200 ) , BACK  < 1200 ) 

REAL  L.AMDA,LNORM(  1200)  ,NOIS(  1200) 

COMMON  FILC0A,F1LC0B,FILC0C,FILC0D,ZETA,VAR,LAMDA,THERMN 
+  ,SH0TN0,G0L.D1  , RXPOW,POWSIG, THRES, RLQAD 
NR- 1200 

K= I NT (  NR/INT  <  G0LD1 ) ) 

NR=IFIX(G0LD1)*K 
XI (0)=0. 

XI (-1 )=0. 

XI < -2 ) -0 . 

Yl  (0)  =  0. 

Y1 <-l)=0. 

Yl ( -2 ) "0. 

*  R1  -  INFORMATION  SIGNAL  (BINARY  1'S  AND  O'S) 

*  XI  ■  GAUSSIAN  DISTRIBUTED  NUMBERS  WHICH  ARE  FED 

*  INTO  THE  LINEAR  FILTER 

*  NOIS  -  GAUSSIAN  DISTIBIJTED  THERMAL  NOISE  CURRENT  VALUES 

*  BACK  =  GAUSSIAN  DISTIBUTED  BACKGROUND  SHOT  NOISE 

*  CURRENT  VALUES 

* 

*  THE  IMSL  ROUTINES  ARE  CALLED  HERE  TO  GENERATE  THESE  VALUES 

CALL  GGUBS(DSEED1 ,NR,R1 ) 

CALL  GGNML.  ( DSEED ,  NR , X 1 ) 

CALL  GGNML ( DSEED2 , NR , NOI S ) 

CALL  GGNML (DSEED3, NR, BACK) 

*  THE  OUTPUT  SEQUENCE  (  Y1(J)  )  FROM  THE  LINEAR  FILTER  IS 

*  DETERMINED  HFRE ,  AS  WELL  AS  THE  VARIANCE  OF  THIS 

*  SEQUENCE  (  AVEVAR  ) 

AVEVAR=0. 

DO  13  J*1,K 

Yl ( J)=ZETA#(FILC0A*X1( J)+3.*FILC0A*X1 ( J-l H 
+3.*FILC0A*X1(  J-  2)fFILC0A*Xl(  J-3 > - ( FILCOC+ 

+3 . *F ILCOD-3 ♦ -F ILCOB ) *Y 1 ( J-l ) - ( 3 . -F ILCOB- 
+FILC0C+3.*FILC0D)*Y1( J-2)-(FILC0B+FILC0D- 
+l.-FILC0C)*Yl(J-3)) 

AVE VAR=AVEVAR+ ( Y 1 ( J ) **2 ) 

13  CONTINUE 

AVEVAR=AVEVAR/FLOAT ( K ) 

DELTAZ* .0001 
VAR1=SQRT (VAR) 

IF ( INDIA.GT . 1 )G0  TO  1 
CALL  DELTA ( VAR , BELT  AY ) 


*  THE  GAUSSIAN  DISTRIBUTED  NUMBERS  ARE  THEN  TRANSFORMED 

*  INTO  LOGNORMALLY  DISTIBUTED  NUMBERS  (  LNORM(J)  )  HERE 

1  DO  61  J=3  ,K 

TARA= ( Y 1 <  J ) /AVEVAR ) -DELTAZ/2 . 

CALL  MDNOR ( T ARA  , T ARA1 ) 

STAR=TARA1*DEL TA2/DELTAY 
CALL  MDNRIS<STAR,STAR1,IER> 

L.NORM  <  J )  =FXP  <  2 .  *VAR1  *  ( STAR1-VAR1 ) )  +DELT AY/2 . 

61  CONTINUE 

*  THIS  SECTION  MULTIPLIES  THE  APPROPRIATE  NUMBER  OF 

*  INFORMATION  BITS  (  R1  <J)  )  BY  THE  PROPER  NUMBER  OF 

*  TURBULENCE  BITS  <  I..MORM <  J )  ),  AND  THEN  ADDS  THE 

*  NOISE  TO  THE  SIGNAL  CURRENT.  THE  RESULTING  CURRENT 

*  LEVEL  IS  THEN  COMPARED  TO  THE  THRESHOLD  LEVEL.  IF 

*  IT  IS  >  THAN  THRESHOLD  A  COUNTER  IS  INTIATED  AND 

*  SlIMED.  THE  TOTAL  SUM  OF  THE  COUTER  VARIABLE  IS  THAN 

*  DETERMINED.  THIS  NUMBER  EQUALS  THE  AMOUNT  OF  l'S 

*  THAT  THE  RECEIVER  'DECIDED'  WERE  SENT. 

NUM~0 

L-1FXX (SOLDI) 

SUM=0. 

DO  1002  J= 1 ,  K 
DO  1003  1  =  1, L 

IF  <  R1 < I +NUM ) . GE . .5)R1 <I+NUM>  =  1 . 

IF(R1 ( I-FNUM) ,LT..5)R1( I+NUM)=0. 

R2  < I +NUM ) =LNORM <  J ) *R1 < I +NUM ) 

NOI S ( I +NUM ) =N0 1 S  < I +NUM ) *SQRT ( THERMN ) 

BACK ( I +NUM ) =BACK <  J  +NUM ) *SQRT  <  SHOTNO ) 
SIGCUR(I+NUM)=SQRT(P0WSIG*R2(I+NUM)/RL0AD> 

TOTAL ( I +NUM ) =SI GCUR ( J  +NUM ) +BACK < I +NUM ) +N0 1 S ( I +NUM ) 
I F  <  TOTAL ( I +NUM ) . GE . THRES )RX= 1 . 

IF  <  TOT  AL  <  J.+NUM  >  .  L  T  .  THRES  )  RX=<> . 

IF(R1 ( I+NUM) . EQ.RX)GO  TO  1003 
SUM=SUM+1. 

1003  CONTINUE 
NUM=NUM+l. 

1002  CONTINUE 
RETURN 
END 

*  SUBROUTINE  'DELTA'  DETERMINES  THE  VALUE  OF  DELTAY  IN 

*  THE  GUJAR/KAVANAUGH  RELATION 

SUBROUT I NE  DELTA ( VAR , BELT AY ) 

DELTAZ= . 0001 
GAUSS=3 » 6B-DELT  AZ/2 . 

VAR1=SQRT (VAR) 

CALL  MDNOR < GAUSS, GAUSS1) 

GAUSS2=GAUSS1*DELTAZ 

DELTAY=0. 

B=EXF'(2  .  *VAR1*  <  3 , 68-VAR1 )  ) 


DO  652  J= If 70000 
Y20=B-DELTAY/2. 

USA= ( 1 . / ( 2 » *VAR1 ) ) #ALOG  <  Y20 ) +UAR1 
CALL  MDNOR<USA,GUN> 
GAUSS4=GAUSS2/GUN 
DIF=GAUSS4-DELTAY 
IF  <  DIF  » l.E » 0  )G0  TO  653 
DELTAY=DELTAY+. 00000015 

652  CONTINUE 

653  CONTINUE 
RETURN 
END 
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